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ON THE INTEGRAL EQUATION OF RENEWAL THEORY 
By Wirity FELLER 


Brown University 


1. Introduction. In this paper we consider the behavior of the solutions of 
the integral equation 


(1.1) u(t) = g(t) + i u(t — x)f(zx) dz, 


where f(t) and g(t) are given non-negative functions.’ This equation appears, 
under different forms, in population theory, the theory of industrial replacement 
and in the general theory of self-renewing aggregates, and a great number of 
papers have been written on the subject.” Unfortunately most of this literature 
is of a heuristic nature so that the precise conditions for the validity of different 
methods or statements are seldom known. This literature is, moreover, abun- 
dant in controversies and different conjectures which are sometimes supported 
or disproved by unnecessarily complicated examples. All this renders an ori- 
entation exceedingly difficult, and it may therefore be of interest to give a 
rigorous presentation of the theory. It will be seen that some of the previously 
announced results need modifications to become correct. 

The existence of a solution u(t) of (1.1) could be deduced directly from a well- 
known result of Paley and Wiener [21] on general integral equations of form 
(1.1).° However, the case of non-negative functions f(¢) and g(t), with which 
we are here concerned, is much too simple to justify the deep methods used by 
Paley and Wiener in the general case. Under the present conditions, the exist- 
ence of a solution can be proved in a simple way using properties of completely 
monotone functions, and this method has also the distinct advantage of showing 
some properties of the solutions, which otherwise would have to be proved 
separately. It will be seen in section 3 that the existence proof becomes most 


natural if equation (1.1)*is slightly generalized. Introducing the summatory 
functions 


a2) vo=[ ua, ro=[ see, Go = [ 2) az, 


1 For the interpretation of the equation cf. section 2. 
* Lotka’s paper [8] contains a bibliography of 74 papers on our subject published before 


be 1939. Yet it is stated that even this list ‘‘is not the result of an exhaustive search.” At 
| the end of the present paper the reader will find a list of 16 papers on (1.1) which have 
| ‘&ppeared during the two years since the publication of Lotka’s paper. 


*This has been remarked also by Hadwiger [3]. 
243 
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equation (1.1) can be rewritten in the form 


(1.3) U(t) = Gd) + [ ‘U(t — 2) dF(2). 


However, (1.3) has a meaning even if F(¢) and G(¢) are not integrals, provided 
F(?) is of bounded total variation and the integral is interpreted as a Stieltjes 
integral. Now for many practical applications (and even for numerical calcula- 
tions) this generalized form of the integral equation seems to be the most 
appropriate one and, as a matter of fact, it has sometimes been used in a more or 
less hidden form (e.g., if all individuals of the parent population are of the same . 
age). Our existence theorem refers to this generalized equation. 
We then turn to one of the main problems of the theory, namely the asymptotic 
behavior of u(t) as t + «. It is generally supposed that the solution u(t) 
“in general’’ either behaves like an exponential function, or that it approaches 


in an oscillating manner a finite limit q; the latter case should arise if f. S() dt = 1, 


thus in particular i in the cases of a stable population and of industrial replace- 

ment. However, special examples have been constructed to show that this is 

not always so.‘ In order to simplify the problem and to get more general condi- 

tions, we shall first (section 4) consider only the question of convergence in mean, 

that is to say, we shall study the asymptotic behavior not of u(t) itself but of 
t 


the mean value u*(t) = ; [ u(x) dx. The question can be solved completely 
o 


using only the simplest Tauberian theorems for Laplace integrals. Of course, 
if u(t) — q then also u*(t) — g, but not conversely. The investigation of the 
precise asymptotic behavior of u(t) is more delicate and ee more refined , 
tools (section 5). 

Most of section 6 is devoted to a study of Lotka’s well-known method of 
expanding u(t) into a series of oscillatory components, and it is hoped that this 
study will help clarify the true nature of this expansion. It will be seen that 
Lotka’s method can be justified (with some necessary modifications). even in 
some cases for which it was not intended, e.g., if the characteristic equation has 
multiple or negative real roots, or if it has only a finite number of roots. On 
the other hand limitations of the method will also become apparent: thus it 
can occur in special cases that a formal application of the method will lead to a 
function u(t) which apparently solves the given equation, whereas in reality it 
is the solution of quite a different equation. 

Of course, most of the difficulties mentioned above arise only when the func- 
tion f(t) has an infinite tail. However, it is known that even computational 
considerations sometimes require the use of such curves, and, as matter of fact, 


‘Cf. Hadwiger [2] and also Hadwiger, ‘‘Zur Berechnung der Erneuerungsfunktion nach 
einer Formel von V. A. Kostitzin,’’ Mitt. Verein. schweizerischer Versich.-Math., Vol. 34 
(1937), pp. 37-43. 
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exponential and Pearsonian curves have been used most frequently in connec- 
tion with (1.1). It will be seen that even in these special cases customary 
methods may lead to incorrect results. Besides, our considerations show how 
much the solution u(t) is influenced by the values of f(t) for £ + ©, and, accord- 
ingly, that extreme caution is needed in practice. The last section contains 
some simple remarks on the practical computation of the solution. 


2. Generalities on equations (1.1) and (1.3). This section contains a few 
remarks on the meaning of our integral equation and on an alternative form 
under which it is encountered in the literature. A reader interested only in the 
abstract theory may pass immediately to section 3. 

Equation (1.1) can be interpreted in various ways; the most important among 
them are the following two: 

(i) In the theory of industrial replacement (as outlined in particular by Lotka), 
it is assumed that each individual dropping out is immediately replaced by a 
new member of zero age. f(t) denotes the density of the probability at the 
moment of installment that an individual will drop out at age t. The function 
g(t) is defined by 


(2.1) a(t) = [ nleyte — 2) ae, 


where (x) represents the age distribution of the population at the moment 
t = 0 (so that the number of individuals of an age between z and x + 6a is 
n(x)dc + o(6x)). Obviously g(t) then represents the rate of dropping out at 
time ¢ of individuals belonging to the parent population. Finally, u(t) denotes 
the rate of dropping out at time ¢ of individuals of the total population. Now 
each individual dropping out at time ¢ belongs either to the parent population, 
or it came to the population by the process of replacement at some moment 
t—2z(0 < x < 2), and hence u(t) satisfies (1.1). It is worthwhile to note that 
in this case 


(2.2) I ” (0) dt = 1, 


since f(t) represents a density of probability. 

(ii) In population theory u(t) measures the rate of female births at time t > 0. 
The function f(t) now represents the reproduction rate of females at age ¢ (that 
is to say, the average number of female descendants born during (t, t + 4t) 
from a female of age t is f(t)dt + o(6t)). If n(x) again stands for the age distri- 
bution of the parent population at t = 0, the function g(t) of (2.1) will obviously 
measure the rate of production of females at time ¢ by members of the parent 
population. Thus we are again led to (1.1), with the difference, however, that 
this time either of the inequalities 


(2.3) | I ” (0) at $1 
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may occur; the value of this integral shows the tendency of increase or decrease 
in the total population. 

Theoretically speaking, f(t) and g(t) are two arbitrary non-negative functions, 
It is true that g(t) is connected with f(t) by (2.1); but, since the age distribution 
n(x) is arbitrary, g(t) can also be considered as an arbitrarily prescribed function. 

It is hardly necessary to interpret the more general equation (1.3) in detail: 
it is the straightforward generalization of (1.1) to the case where the increase or 
decrease of the population is not necessarily a continuous process. This form 
of the equation is frequently better adapted to practical needs. Indeed, the 
functions f(t) and g(t) are usually determined from observations, so that only 
their mean values over some time units (years) are known. In such cases it is 
sometimes simpler to treat f(t) and g(t) as discontinuous functions, using 
equation (1.3) instead of (1.1). For some advantages of such a procedure 
see section 7. It may also be mentioned that the most frequently (if not the 
only) special case of (1.1) studied is that where g(t) = f(t). Now it is apparent 
from (2.1) that this means that all members of the parent population are of 
zero age: in this case, however, there is no continuous age-distribution 7(z). 
Instead we have to use a discontinuous function n(x) and write (2.1) in the form 
of a Stieltjes integral. Thus discontinuous functions and Stieltjes integrals 
present themselves automatically, though in a somewhat disguised form, even 
in the simplest cases. 

At this point a remark may be inserted which will prove useful for a better 
understanding later on (section 6). In the current literature we are frequently 
confronted not with (1.1) but with 


(2.4) ui) = [ ” wt ~ adds, 


together with the explanation that it is asked to find a solution of (2.4) which 
reduces, for ¢ < 0, to a prescribed function A(t): Now such a function, as is 
known, exists only under very exceptional conditions, and (2.4) is by no means 
equivalent to (1.1). The current argument can be boiled down to the following. 
Suppose first that the function g(t) of (1.1) is given in the special form 


(2.5) g(t).= f° re — xt) ae, 


where h(x) is a non-negative function defined for x < 0. Since the solution 
u(t) of (1.1) has a meaning only for t > 0, we are free. to define that u(—t) = 
h(—t) fort > 0. This arbitrary definition, then, formally reduces (1.1) to (2.4). 
It should be noted, however, that this function u(#) does not, in general, satisfy 
(2.4) for t < 0, for h(t) was prescribed arbitrarily. Thus we are not, after all, 
concerned with (2.4) but with (1.1), which form of the equation is, by the way, 
the more general one for our purposes. If there really existed a solution of 
(2.4) which reduced to A(t) for t < 0, we could of course define g(¢) by (2.5) and 
transform (2.4) into (1.1) by splitting the interval (0, ©) into the subintervals 
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(0, é) and (¢, ©). However, as was already mentioned, a solution of the required 
kind does not exist in general. It will also be seen (section 6) that the true 
nature of the different methods and the limits of their applicability can be under- 
stood only when the considerations are based on the proper equation (1.1) and 
not on (2.4). 


3. Existence of solutions. 
THEOREM 1. Let F(t) and G(t) be two finite non-decreasing functions which 
are continuous to the right’. Suppose that 


(3.1) F(0) = GO) = 0, 
and that the Laplace integrals® 


(3.2) o(s) = [ “edF(), —-y(s) = [ ” aan 


converge at least fors >oa>0'. Incase that lim ¢(s) > 1, let o’ > a be the root* 
s—o+0 
of the characteristic equation ¢(s) = 1; in case lim ¢(s) < 1, puto’ =o. 
s—ot+0 
Under these conditions there exists for t > 0 one and only one finite non-decreasing 


function U(t) satisfying (1.3). With this function the Laplace integral 


(3.3) én f * 6 dU) 


5 It is needless to emphasize that this restriction is imposed only to avoid trivial am- 
biguities. 
* The integrals (3.2) should be interpreted as Lebesgue-Stieltjes integrals over open 
intervals; thus 
oo 
¢g(s) = lim e~* dF(t), 
e—-+0 J, 
which implies that ¢(s) ~ 0 ass— «. Alternatively it can be supposed that F(t) and 
G(t) have no discontinuities at ¢ = 0. Continuity of F(t) at t = 0 means that there is no 
reproduction at zero age. This assumption is most natural for our problem, but is by no 
means necessary. In order to investigate the case where F(t) has a saltusc > 0 att = 0, 
one should take the integrals (3.2) over the closed set [0, ~], so that 
oO 
g(s) = c+ lim e~* dF(t). 
e—+0 


It is readily seen that Theorem 1 and its proof remain valid if 0 < ¢c <1. However, if 
c > 1, then (1.3) plainly has no solution U(t). The continuity of G(t) at ¢ = 0 is of no 
importance and is not used in the sequel. 

7 The condition is formulated in this general way in view of later applications (cf., e.g., 
the lemma of section 4). In all cases of practical interest o = 0. 

® »(s) is, of course, monotonic for s > o and tends to zeroass—> ~. In order to ensure 
the existence of a root of ¢(s) = 1, it is sufficient to suppose that the saltus c of F(t) at t = 0 
is less than 1 (cf. footnote 6). 
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converges for s > a’, and 


(3.4) v(s) = 7. 


Proor: A trivial computation shows that for any finite non-decreasing solu- 
tion U(t) of (1.3) and any T > 0 we have 








T - * 7 T2 
—st U _ —st —62 —st ™ 
I e* due) I e ac) + | é ar (2) | e* dU(0); 
herein all terms are non-negative and hence by (3.2) 
7 , 
[ eau < v0) + 06) [ eave. 
0 


Now ¢(s) < 1 for s > o’, and hence it is seen that the integral (3.3) exists for 
s > o’ and satisfies (3.4). On the other hand it is well-known that the values of 
w(s) for s > o’ determine the corresponding function U(t) uniquely, except for 
an additive constant, at all points of continuity. However, from (1.3) and (3.1) 
it follows that U(0) = 0 and, since by (1.3) U(¢) is continuous to the right, the 
monotone solution U(t) of (1.3), if it exists, is determined uniquely. 

To prove the existence of U(¢) consider a function w(s) defined for s > a’ by 
(3.4). It is clear from (3.2) that ¢(s) and y(s) are completely monotone func- 
tions, that is to say that ¢(s) and y(s) have, for s > o, derivatives of all orders 
and that (—1)"9‘”(s) > 0 and (—1)"y‘"(s) > 0. Wecan therefore differentiate 
(3.4) any number of times, and it is seen that w”’(s) is continuous for s > o’. 
Now a simple inductive argument shows that (—1)"w'”(s) is a product of 
{1 — y(s)}~"*” by a finite number of completely monotone functions. It 
follows that (—1)"w”(s) > 0, so that w(s) is a completely monotone function, 
at least for s > o’. Hence it follows from a well-known theorem of S. Bernstein 
and D. V. Widder’ that there exists a non-decreasing function U(t) such that 
(3.3) holds for s > o’. Moreover, this function can obviously be so defined that 
U(0) = 0 and that it is continuous to the right. Using U(t) let us form a new 
function 


(3.5) Vi) = [ ' U(t — 2) dF(2). 


V(t) is clearly non-negative and non-decreasing. It is readily verified (and, of 
course, well-known) that 


vs) = [ e*aV® = w()o(. 
0 
It follows, therefore, from (3.4) that ¥(s) = w(s) — y(s), and this implies, by the 


* This theorem has been repeatedly proved by several authors; for a recent proof cf. 
Feller [19]. 
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uniqueness theorem for Laplace transforms, that V(t) = U(#) — G(t). Combin- 
ing this result with (3.5) it is seen that U(¢) is a solution of (1.3). 

THEOREM 2. Suppose that f(t) and g(t) are measurable, non-negative and 
bounded in every finite intervralO0Q < t < T. Let the integrals 


3.6) ols) = [ “etd, (8) = [ ” g(t) dt 


converge fors > a. Ther there exists one and only one non-negative solution u(t) 
of (1.1) which is bounded in every finite interval”. With this function the integral 


(3.7) o(@) = f ” u(t) dt 


converges at least for s > o’, where o’ = o if lim ¢(s) < 1, and otherwise a’ > @ 


is defined as the root of the characteristic equation y(s) = 1. Fors > o' equation 
(3.4) holds. 

If f(t) is continuous except, perhaps, at a finite number of points then u(t) — g(t) 
is continuous. 

Proor: Define F(¢#) and G(¢) by (1.2). Under the present conditions these 
functions satisfy the conditions of Theorem 1, and hence (1.3) has a non-decreas- 
ing solution U(t). Consider, then, an arbitrary interval 0 < ¢ < T and suppose 
that in this interval f(t) < M and g(t) << M. If0 <i<t+h< T we have 
by (1.3) 


o<i{ue+aA) — UW} 
= 1 (Go +h) — G@} + if” Ult + h — 2)f(2) de 
h h J: 


+} [ (uG+h—2) — Ub —2)\fle) dz 
0 
<m+mur) + [ (ue+h—2) - U(t — x)} dz 


- tth 
=u mun) +M [veya —™ [ varey 
<M + 2MU(T). 


Thus U(t) has bounded difference ratios and is therefore an integral. The 
derivative U’(t) exists for almost all tand 0 < U’(t) < M. Accordingly we can 
differentiate (1.3) formally, and since U(0) = 0 it follows that u(é) = U’(t) 
satisfies (1.1) for almost all ¢. However, changing u(¢) on a set of measure zero 
does not affect the integral in (1.1), and since g(¢) is defined for all ¢ it is seen that 


10 Without the assumptions of positiveness and boundedness this theorem reduces to a 
special case of a theorem by Paley and Wiener (21]; cf. section 1, p. 243. 
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u(t) can be defined, in a unique way, so as to satisfy (1.1) and obtain (1.3), 
Since the solution of (1.3) was uniquely determined it follows that the solution 
u(t) is also unique. Obviously equations (3.7) and (3.3) define the same function 
w(s), so that (3.4) holds, and (3.7) converges for s > o’. 

Finally, if f(t) has only a finite number of jumps, the continuity of u(t) — g(t) 
becomes evident upon writing (1.1) in the form 


t 
u(t) — g(t) = I u(x)f(t — x) dz. 
4. Asymptotic properties. In this section we shall be concerned with the 
asymptotic behavior as t — © not of u(t) itself but of the mean value u*(t) = 
t 
; [ u(r)dr. If u(t) tends to the (not necessarily finite) limit C, then obviously 
() 


also u*(t).— C, whereas the converse is not necessarily true. For the proof of 
the theorem we shall need the following obvious but useful 
Lemma: If u(t) > 01s a solution of (1.1) and of 


(4.1) u(t) =eut), AQ =&fO, git) = gd, 
then u;(t) is a solution of 


ud) = ad) + f “walt — 2)fa) de. 


THEOREM 3: Suppose that using the functions defined in Theorem 2 the integrals 


(4.2) I f(t) dt = a, I g(t) dt = b, 
are finite. 
(i) In order that 
. sae 1 , 
(4.3) ur(t) = 2 I uls)de + C 


as t— ©, where C is a positive constant, it is necessary and sufficient that a = 1, 
and that the moment, 


(4.4) I  (iet<w 


be finite. In this case 





(4.5) C= 
(ii) If a < 1 we have 


. b 
(4.6) [ a= —. 


(1.8), 
lution 
nction 


— g(t) 
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(iii) If a > 1 let o’ be the positive root of the characteristic equation y(s) = 1 
(cf. (3.2)) and put 


" "tt £(t) dt = m. 
(4.7) [ e f@® ™m 
Then 
i a b 
(4.8) lim— | e° u(r) dr = —. 
t 0 m 


t-*0 

RemaRK: The case a = 1 corresponds in demography te a population of 
stationary size. In the theory of industrial replacement only the case a = 1 
occurs; the moment m is the average lifetime of an individual. The case a > 1 
corresponds in demography to a population in which the fertility is greater than 
the mortality. As is seen from (4.8), in this case the mean value of u(é) increases 
exponentially. It is of special interest to note that in a population with a < 1 
the integral (4.6) always converges. 

Proor: By (4.2) and (3.7) 


(4.9) lim ¢(s) = a, lim y(s) = b. 
s—+0 s—+0 


If a < 1, it follows from (3.4) that lim w(s) = b/(1 — a) is finite. Since u(t) >0 
this obviously implies that (4.6) holds. This proves (ii). 
If a = 1 and mis finite, it is readily seen that 


lim 1 — os) o(s) = 


m, 
s—+0 8 


and hence by (3.4) 
ins) = i, 0) in, y= 
By a well-known Tauberian theorem for Laplace integrals of non-negative 
functions” it follows that u*(t) > =. Conversely, if (4.3) holds it is readily 
seen that” 
11 (4.2) implies the finiteness of m, . 


2 Cf. e.g. Doetsch [18], p. 208 or 210. 


18 Indeed, if (4.3) holds and if U(t) is defined by (1.2), then there isa M = M(e) sch that 
| Ut) —Ct|<M+e. Now 


eo 
¢(s) = | e~*U(t) dt, 
0 
and hence 


ove) -C a” I e-*(U() — C2) dt, 
0 


| se(s) — C] <#/ e"(M + d) dt =sM +.«. 
0 
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lim sw(s) = C, 
s—-+0 

which in turn implies by (3.4) and (4.9) that 


.- 1—g(s)_ b 
-~ 8 - c 
This obviously means that the moment (4.4) exists and equals b/C. This: 
proves (i). 
Finally case (iii) reduces immediately to (ii) using the above lemma with 
k = —o’. This finishes the proof. 
It may be remarked that the finiteness of the integrals (4.2) is by no means. 
necessary for (4.3). This is shown by the following 
EXAMPLE: Let 





i) = me g(t) = em 
It is readily seen that with these functions a = 1, butb = ©. Now” ¢(s) = 
e~V* and y(s) = ¢ V*/+/s, so that 
eve 
Vs(1 —ev*)’ 


Thus sw(s) — 1 as s —> + 0, and hence u*(t) 1. In this particular case it can 
even be shown that the solution u(t) itself tends to 1 ast — ©. 

In practice, however, the integrals (4.2) will always exist, and accordingly we 
restrict the consideration to this case. 


w(s) = 


5. Closer study of asymptotic properties. In this section we shall deal almost 
exclusively with the most important special case, namely where 


(5.1) [ soa =1. 

0 
The question has been much discussed whether in this case necessarily u(t) — C 
as t — ©, which statement, if true, would be a refinement of (4.3). Hadwiger 
[2] has constructed a rather complicated example to show that u(t) does not 


necessarily approach a limit. Now this can also be seen directly and without 
any computations. Indeed, if u(¢) — C and if (5.1) holds, then obviously 


t 
lim [ u(t — 2)f(2) dr = ©, 
t-re 40 
and hence it follows from (1.1) that g(t) +0. In order that u(t) —C it is therefore 


14 The integrals can be evaluated by elementary methods, and are known; ef. e.g. 
Doetsch [18], p. 25. 
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necessary that g(t) —> 0, and this proves the assertion. In Hadwiger’s example 
lim sup g(t) = ©, which makes his computations unnecessary. 

It can be shown in a similar manner that not even the condition g(t) — 0 is 
sufficient to ensure that u(t) ~ C. Some restriction as to the total variation of 
f(t) seems both necessary and natural (conditions on the existence of derivatives 
are not sufficient). In the following theorem we shall prove the convergence of 


This: u(t) under a condition which is, though not strictly necessary, sufficiently wide 
to cover all cases of any possible practical interest. 

vith THEOREM 4: Suppose that with the functions f(t) and g(t) of Theorem 2 

saa (5.2) I fit) dt =1, I g(t) dt =b < . 


Suppose moreover that there exists an integer n = 2 such that the moments 
(5.3) m, = [ #*#(t) dt, beh Si +~.% 
o 


) = are finite, and that the functions f(t), tf(t), Uf(t), --- , t” °f(t) are of bounded total 
variation over (0, ©). Suppose finally that 


(5.4) lim t"*g()=0 and limt™ [ g(z)dx =0. 
t—*0 t—*c t 
Then 
can 
(5.5) lim u(t) = © 
we tree = 
and 
a b \ 
10st (5.6) lim t” “4 u(t) — —? = 0. 
tc ™m j 
Remark: As it was shown in section 4, the case where | f(idt>1 
0 
can readily be reduced to the above theorem by applying the lemma of section 4 
oC with k = o’, where a’ is the positive root of ¢(s) = 1: it is only necessary to 
| suppose that e~* ‘f(t) is of bounded total variation and that e~* ‘g(t) +0. Ob- 
= viously all moments of e~” ‘f(t) exist, so that the above theorem shows that 
a u(t) = e~*' ‘u(t) tends to the finite limit b’/m; , where 
b! = I e*'' g(t) dt, m, = | e* ‘if (t) dt. 
0 
Thus in this case and under the above assumptions u(t) ~ oat e”'' so that the 
1 
ore renewal function increases exponentially as could be expected. If however 


[1 dt <1, 





254 WILLY FELLER 


u(t) will in general not show an- exponential character. If f(t) is of bounded 
variation and has a finite moment of second order, and if g(t) — 0, then it can be 
shown that u(t) —- 0. However, the lemma of section 4 can be applied only if 
the integral defining ¢(s) converges in some negative s-interval containing a value 
s’ such that ¢(s’) = 1, and this is in general not the case. 

Proor: The proof of Theorem 4 will be based on a Tauberian theorem due to 
Haar’. With some specializations and obvious changes this theorem can be 
formulated as follows. . 

Suppose that /(t) is, for ¢ > 0, non-negative and continuous, and that the 
Laplace integral 


(5.7) id = [ ” &*1(t) dt 


converges for s > 0. Consider \(s) as a function of the complex variables = 
x + iy and suppose that the following conditions are fulfilled: 

(i) For y ¥ 0 the function \(s) (which is always regular for x > 0) has con- 
tinuous boundary values A(iy) as rx — +0, for z > 0 andy #0 


(5.8) na) = © + vie), 


where (iy) has finite derivatives y’(iy), --- , vw (iy) and y"(iy) is bounded 
in every finite interval; 


(ii) [. Na + ty) dy 


converges for some fixed xz > 0 uniformly with respect tot > T > 0; 
(iii) A(x + ity) ~0 as y > +, uniformly with respect to x > 0; 
(iv) d’(iy), (iy), --- ; A (dy) tend to zero as y > +; 

(v) The integrals 


ae ae 
[ e*) (iy) dy and / (iy) dy 
Cs) v2 
(where y; < 0 and y2 > 0 are fixed) converge uniformly with respect tot > T > 0. 
Under these conditions 
(5.9) lim t{l(t) — C} = 0. 


Now the hypotheses of this theorem are too restrictive to be applied to the 
solution u(t) of (1.1). We shall therefore replace (1.1) by the more special 
equation 





(5.10) v(t) = A(t) + I '~ ae. 


16 Haar [20] or Doetsch [18], p. 269. 
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where 
(5.11) h(t) = [x0 — 2x) f(x) dz. 


Plainly Theorem 2 can be applied to (5.10). It is also plain that h(t) is bounded 
and non-negative and that (by (5.1)) 


(5.12) [ h(t) dt = 1, 


(5.13) x(s) = I e*'h(t) dt = ¢(s). 
Accordingly we have by Theorem 2 


st ” —st i ¢'(s) 
(5.14) ¢(s) = I ev) dt = or. 


We shall first verify that ¢(s) satisfies the conditions of Haar’s theorem with 
r=n-— 2. For this purpose we write 
(5.15) SO = AO — fh, 
where fi(¢) and f2(¢) are non-decreasing and non-negative functions which are, 
by assumption, bounded: 
(5.16) 0<fif) < M, 0< f(t) < M. 


(a) We show that v(t) is continuous. Now by Theorem 2 the solution v(é) 
of (5.10) is certainly continuous if h(é) is continuous; however, that h(t) is con- 
tinuous follows directly from (5.11) and the fact that the functions 


[ne -aseyar ond [ple aye) az 


are continuous. 
(b) In view of (5.1) the function 9(s) exists for z = R(s) > 0. Obviously 
|e(z +iy)|<1forz>0. Now 


1- oy) = [ae poae 


= [ (1 — cos yt) f(t) dt +7 [ sin yt-f(t) dt, 


and, since 1 — cos yt > 0 and f(t) > 0, the equality g(iy) = 1 for y + 0 would 
imply that f(t) = 0 except on a set of measure zero. It is therefore seen that 
g(x + ty) ¥ 1 forall z > O and forz = 0, y ¥ 0. 

It follows furthermore from (5.3) that fork = 1, --- ,n and xz > 0 the deriva- 
tives 


e”(s) = [ (— t)* e** f(t) at 








256 WILLY FELLER 
exist and that 
jim eo” (a + ty) = o (iy). 
Finally, it is readily seen that in the neighborhood of y = 0 we have 
e(iy) = [po at 


(5.17) =1—miy+ 5 (iy) metas 





+ (—1)"" ae (iy)"* + O(ly |"). 


(c) From what was said under (b) it follows by (5.14) that ¢(s) is regular for 
xz > 0, and that ¢(s), ¢’(s), --- , ¢°(s) approach continuous boundary values 
‘ as § = x + ty approaches a point of the imaginary axis other than the origin. 
Now put 


(5.18) jio~ 28 2, 





so that by (5.14) 
(5.19) i «+960. 
mS 


For z > 0 and z = 0, y ¥ O the function ¥(z + iy) is obviously continuous; 
the derivatives (iy), --- ,v'” (dy) -exist. To investigate the behavior of 
¥(zy) in the neighborhood of y = 0 put 


(5.20) Ply) = m — Fy) + — - — (1 i. 





By (5.17), (5.18) and (5.20) 


_\ | {1 — iyP(y)}? 1141 et 

(5.21) vy) = | PaO 112 + ocyir. 
Now the expression in brackets represents an analytic function of y which 
vanishes at y = 0. Hence y¥(iy) = B(y) + O(|y|"”), where B(y) denotes a 
power series. It follows that the derivatives y/(iy), --- , y°"-? (iy) exist for all 
real y (including y = 0) and are bounded for sufficiently small | y | : since they 
are continuous functions they are bounded in every finite interval. 

(d). Next we show that there exists a constant A > 0 such that for sufficiently 
large | y | . 


(5.22) lol + | < 


uniformly inz >0. By (5.15) 


(5.23) ¢(s) = I ° {cos yt — isin yt}e *‘{ fi(t) — fo(t)} dt. 
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Now f(t) is non-decreasing and accordingly by the second mean-value theorem 
we have for any T > O and y 


? r ; ; 
I cos yt-fi(t) dt = fi(T) I cos yedt = #(7) Y= STH, 


where 7 is some value between 0 and 7' (depending, of course, on y; at points of 
discontinuity, f:(7’) should be replaced by lim f,(t)). Hence by (5.16) 
t—T—0 


| cos yt-e ~-f,(t) ar < 2M 
0 ly | 


Treating the other terms in (5.23) in a like manner, (5.22) follows. 
Combining (5.22) with (5.14) it is seen that for sufficiently large | y | 

2a" 

y° 

uniformly in z > 0. This shows that the assumptions (ii) and (iii) of Haar’s 
theorem are satisfied for \(s) = ¢(s). In order to prove that also conditions 
(iv) and (v) are satisfied it suffices to notice that the proof of (5.22) used only 
the fact that f(t) is of bounded total variation. Now ¢“’(s) is the Laplace trans- 


form of (—t)*f(t), and, since ¢‘f(t) is of bounded total variation for k < n — 2, 
it follows that 


|(s)| < 


|~(s) | = Oly |”), k= 1, 2,---,n — 2, 


for sufficiently large | y |, uniformly inz > 0. Differentiating (5.14) k times it is 
also seen that 


Is(s) | = Oy I>), k= 1,2,---,n —_ 2, 
asy— + ©, uniformly with respect to z > 0. 
This enumeration shows that v(s) = l(t) and A(s) = ¢(s) satisfy all hypotheses 
of Haar’s theorem with r = n — 2andC = 1/m,. Hence 
(5.24) lim {oto ~ 1} = 0. 
t—*co m™m 
Returning now to (5.14) we get 
w(s) = y(s) + y(s)e(s) + r(s)f(s), 
or, by the uniqueness property of Laplace integrals, 


(5.25) u(t) = g(t) + [ g(x) f(t — x) dx + [ g(x)o(t — x) dz 


= g(t) + u(t) + w(t) 


(which relation can also be checked directly using (5.10)). Let us begin with 
the last term. We have by (5.2) 
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u(t) — - = [ g(t — 2x) {ote - 1 hae, 
and hence 


eo 


t 
u(t) — | < 2"° / g(t — x)x"” dx 
m t/2 











v(x) — 2. 


t 
n—2 
+t | . g(y) 





1| 
v(t — y) on, |: 


If t is sufficiently large we have by (5.24) in the first integral x” ” 





ua) 1] <e 


In the second integral v(t — y) — -. is bounded, and hence by (5.4) 
1 


. —2 
lim ¢” 
t-—*c 





b 
u(t) — >| = 0. 

m 
The same argument applies (even with some simplifications) also to the second 
term in (5.24); it follows that 


lim t” u(t) = 0, 


whilst t” *g(t) + 0 by assumption (5.4). Now the assertion (5.6) of our theorem 
follows in view of (5.25) if the last three relationships are added. This finishes 
the proof of Theorem 4. 

It seems that the solution u(t) is generally supposed to oscillate around its 
limit b/m, ast — ©. It goes without saying that such a behavior is a priori 
more likely than a monotone character. It should, however, be noticed that 
there is no reason whatsoever to suppose that u(t) always oscillates around its 
limit. Again no computation is necessary to see this, as shown by the following 

EXxaMPLeE: Differentiating (1.1) formally we get 





u(t) = g(t) + gO) f@ + [ u(t — x) f(x) dz, 


which shows that, if g(t) and f(t) are sufficiently regular, w’(t) satisfies an integral 
equation of the same type as u(t). Thus if 


g(t) + gO)f) 2 0 


for all ¢, we shall have u’(t) > 0, and u(t) is a monotone function. In particular, 
if g’(t) + g(O)f() = 0, then u’(¢t) = 0 and u(t) = const. For example, let f(t) = 
g(t) =e‘. Then g(s) = y(s) = 1/(s + 1) and hence w(s) = 1/s, which is the 
Laplace transform of u(t) = 1. It is also seen directly that u(t) = 1 is the 
solution. We have however the following 

Turorem 5°: If the functions f(t) and g(t) of Theorem 4 vanish identically for 
t > T > 0, then the solution u(t) of (1.1) oscillates around its limit b/mast— ~. 


16 Under some slight additional hypotheses and with quite different methods this theorem 
was proved by Richter [16]. 
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Proor: For ¢t > T equation (1.1) reduces to 
t 
u(t) = [u(t — 2)f(@) de, 
t—T 


t 
and since | f(x) dx = 1 it follows that the maxima of u(¢) in the intervals nT < 
t—T 


t < (n + 1)T form, for sufficiently large integers n, a non-increasing sequence. 
Similarly the corresponding minima do not decrease. Since u(t) — b/m,, by 
Theorem 4, it follows that the minima do not exceed b/m, and the maxima are 
not smaller than b/m . 


6. On Lotka’s method. Probably the most widely used method for treating 
equation (1.1) in connection with problems of the renewal theory is Lotka’s 
method. As a matter of fact this method consists of two independent parts. 
The first step aims at obtaining the exact solution of (1.1) in the form of a series 
of exponential terms (this is achieved by an adaptation of a method which was 
used by P. Herz and Herglotz for other purposes. The second part of Lotka’s 
theory consists of devices for a convenient approximative computation of the 
first few terms of the series. While restricting ourselves formally to Lotka’s 
theory, it will be seen that some of the following remarks apply equally to other 
methods. 


Lotka’s method rests essentially on the fundamental assumption that the 
characteristic equation 


(6.1) o(s) = 1 


has infinitely many distinct simple” roots s , s:, --- , and that the solution u(t) 
of (1.1) can be expanded into a series 


(6.2) u(t) = x A,e** 


where the A; are complex constants. The argument usually rests on an assumed 
completeness-property of the roots. Thus, starting from (2.4) it is required that 
(6.2) reduces to h(t) for ¢ < 0; in other words, that an arbitrarily prescribed 
function h(x) be, for z < 0, respresentable in the form 


(6.3) h(x) = X A,e’** (x <0). 


In practice we are, of course, usually not concerned with h(t) but with g(t) (cf. 
(2.5)), and according to Lotka’s theory the coefficients A; of the solution (6.2) 
of (1.1) can be computed directly from g(t) in a way similar to the computation 
of the Fourier coefficients. 


Lotka’s method is known to lead to correct results in many cases and also to 





1” Hadwiger [3] objected to the assumption that all roots of (6.1) be simple. The modifi- 
cations which are necessary to cover the case of multiple roots also will be indicated below. 
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have distinct computational merits. On the other hand it seems to require a 
safer justification, since its fundamental assumptions are rarely realized. Thus 
clearly an arbitrary function A(x) cannot be represented in the form (6.3): to 
see this it suffices to note that (6.1) frequently has only a finite number of roots 
(cf. also below). It should also be noted that, the series (6.3) having regularity 
properties as are assumed in Lotka’s theory, any function representable in the 
form (6.3) is necessarily a solution of the integral equation (2.4), whereas the 
theory requires us to construct a solution u(t) which reduces to an arbitrarily 
prescribed function h(t) for ¢ < 0, (which frequently is an empirical function, 
determined by observations). Nevertheless, it is possible to give sound founda- 
tions to Lotka’s method so that it can be used (with some essential limitations 
and modifications) sometimes even in cases for which it originally was not 
intended. For this purpose it turns out to be necessary that all considerations 
be based on the more general equation (1.1), instead of (2.4) (ef. also section 2). 

Before proceeding it is necessary to make clear what is really meant by a root of 
(6.1). The function ¢g(s) is defined by (3.2), and the integral will in general 
converge only for s-values situated in the half-plane R(s) > o. Usually only 
roots situated in this half-plane are considered. It is also argued that ¢(s) 
is, for real s, a monotone function, so that (6.1) has at most one real root: ac- 
cordingly the terms of (6.2) are called “oscillatory components.’’ However, 
the function ¢(s) can usually be defined by analytic continuation even outside 
the half-plane R(s) > o, and, if this is done, (6.1) will in general also have roots 
in the half-plane R(s) <o. It will be seen in the sequel that these roots play 
exactly the same role for the solution u(t) as the other ones, and that the ap- 
plicability of Lotka’s method depends on the behavior of ¢(s) in the entire 
complex s-plane. It may be of interest to quote an example where (6.1) has 
infinitely many real and no other roots. 

Example”: Let 


1 1 
(6.4) f®) = 2a pr° , t>0; 


18 This was stated in particular by Hadwiger [3] and Hadwiger and Ruchti [6]; accord- 
ingly the results of the latter paper (obtained by methods quite different from Lotka’s) 
need some modifications. 

19Cf. the example at the end of section 4. A function closely related to (6.4) 
plays an important role in two recent papers by Hadwiger [4] and [5]. Hadwiger’s conclu- 
sion, if it could be justified, would fundamentally change the aspect of the whole theory. 
The conclusion reached by Hadwiger seems to be that for any biological population the 
reproduction function should be of the form u(t) = Zu,(t), where un(t) represents the 
contribution of the nth generation and 


(*) uy(t) = ee e~At+Can—n2a2/t 


Vr {3/2 


Here a, A andC areconstants. Clearly (*) is a generalization of (6.4). Now his conclusion 
is based on the arbitrary assumption that u,(t) should be of the form u(t) = ¥(z, na) 
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It is easily seen that g(s) = e~V*. The integral (3.2) converges only for R(s) > 
0, but ¢(s) is defined as a two-valued function in the entire s-plane. The roots of 
(6.1) are obviously 5, = —4 k’x’, so that all of them are real and simple. If 
g(t) = f(t), we get by (3.4) 


eve = 
= —nVe 
w(s) ooo de ' s real, > 0. 


Now e-"v* is the Laplace transform of —~— ¢~"*/ 





, and hence it is readily 





24/x #2 
seen that the solution u(¢) can be written in the form 
_ 1 = —n2/4t , 
(6.5) u(t) — 2./m {3/2 x ne ’ 


of course, this expansion is not of form (6.2) and shows no oscillatory character. 

From now on we shall consistently denote by ¢(s) the function defined by the 
integral (3.4) and by the usual process of analytic continuation; accordingly we 
shall take into consideration all roots of (6.1). The main limitation of Lotka’s 
theory can then be formulated in the following way: Lotka’s method depends 
only on the function g(¢) and on the roots of (6.1). Now two different functions 
f(t) can lead to characteristic equations having the same roots. Lotka’s method 
would be applicable to both onl:’ if the corresponding two integral equations 
(1.1) had the same solution u(t). This, however, is not necessarily the case. 
Thus, if Lotka’s method is applied, and if all computations are correctly per- 
formed, and if the resulting series for u(t) converges uniformly, there is no 
possibility of telling which equation is really satisfied by the resulting u(é): 
it can happen that one has unwittingly solved some unknown equation of type 
(1.1) which, by chance, leads to a characteristic equation having the same roots 
as the characteristic equation of the integral equation with which one was really 
concerned. Indeed this happens in the following example which is familiar in 
connection with our problem. It is illustrative also for other purposes: thus it 
shows not only limitations of Lotka’s method, but also that this method can be 
modified so as to become applicable in some cases where the characteristic equa- 
tion has only a finite number of roots. 


where ¥(z, a) is independent of n. To my mind Hadwiger’s result shows only the im- 
practibility of this axiom. However, Hadwiger’s result is not correct even under his assump- 
tion. Indeed, he derives for y(z, a) the functional equation 


t 
(+) v(z, a+b) = I ve — & av(E, b) dé, 
0 


which is well-known from the theory of stochastic processes. Now Hadwiger merely 
verifies the known result that (*) leads to a solution of (**). However, (**) has infinitely 
many other solutions (it is possible to write down expressions for their Laplace transforms, 
although it is difficult to express the solutions themselves explicitly). This, of course, 
renders Hadwiger’s result illusory. 
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ExamPLe: Pearson type III-curves.° Consider the integral equation (1.1) 
in the following two cases: 


(D) 0) = od = A) = aq tite 
and 
(IT) f(t) = g() = fult) = 4 e™*. | 
It is readily seen (and well known) that the corresponding Laplace transforms are 
(1) gr(s) = Ep 
and 
(ID) gu(s) = ail 
(s + 1)° 


respectively. Thus in both cases the characteristic equation has the same roots, 
namely 


3 pn 
s, = 0, 82,3 = “3 & 5 V3; 


of which only the first one lies in the half-plane of convergence of the integral 
(3.4). Lotka’s method is not applicable since there are only three roots. How- 
ever, in the second case, an expansion of type (6.2) is possible. Indeed, we have 
by (3.4) 

ou(s)  _ 1 
1—gu(s) s*+3s?+ 3s 

; —_™ : 5 + 

a 6 23 i avi 3 
Sete bVv3 et 5+ iA 


wy(s) = 





—at 


now 1/(s + a) is the Laplace transform of e 
u(t) in the form 


, and hence we obtain the solution 





1 1 a 4(—8+iV3)¢ ( a ) 4(—3-iV3) t 
u(t) = - —|- — —=le —{-+ —je 
uO = 3 (; =) 6 2/3 
a 1 pati cos v3, _ i e? sin mee 


20 General Pearson curves have been investigated recently in connection with (1.1) by 
Brown [1], Hadwiger and Ruchti [6] and Rhodes [15]. Hadwiger and Ruchti use a method 
of their own, but they are also led to the study of the characteristic equation (6.1) in 4 
slightly disguised form: their result needs a modification since they arbitrarily drop the 
roots lying in the halfplane of divergence of the integral ¢(s). 


(1.1) 


ms are 


itegral 
How- 
e have 


»lution 


(1.1) by 
method 
1) in a 
rop the 
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which is an expansion of type (6.2). In the first of the above examples we get 
for real positive s 


gi(s) =< 1 
1 — ¢i(s) = (s + 1)%" 
and it is readily seen that this is the Laplace transform of the solution 


wy(s) = 


— = 1 3(n—2) /2 
u(t) nm) I'(3n/2) . 


The series is convergent for ¢ > 0, but obviously this solution cannot be repre- 
sented in a form similar to (6.2). 


A similar remark applies to the general Pearson-type III curve 


f@® = Ate, 
where A, a, 8 are positive constants; the corresponding Laplace transform is 
1 
¢(s) = AT(8 + 1) (¢+ayPH" 


These preparatory remarks enable us to formulate rigorous conditions for the 
existence of an expansion of type (6.2). The following theorem shows the limits 
of Lotka’s method, but at the same time it also represents an extension of it. 
In the formulation of the theorem we have considered only the case of absolute 
convergence of (6.2). This was done to avoid complications lacking any practi- 
cal significance whatsoever. The conditions can, of course, be relaxed along 
customary lines. 

THEOREM 6: In order that the solution u(t) of Theorem 2 be representable in 
form (6.2), where the series converges absolutely for t > 0 and where the s, denote the 
roots of the characteristic equation” (6.1), it is necessary and sufficient that the La- 
place transform w({s) admit an expansion 


(6.6) dhe MO «ES 


1 — ¢(s) 8 — & 
and that =| A; | converges absolutely. The coefficients A; are determined by 








(Sx) 
6.7 Ayo ~ Le. 
(6.7) . (Sk) 
In particular, it is necessary that w(s) be a one-valued function.” 

Proor: All roots s of (6.1) satisfy the inequality R(s.) < o’, where o’ was 
defined in Theorem 2. It is therefore readily seen that in case =| A; | con- 
verges, the Laplace transform of (6.2) can be computed for sufficiently large 


*1 The number of roots may be finite or infinite. It should also be noted that it is not 
required that s, > ©. If the s, have a point of accumulation, w(s) will have an essential 
singularity. That this actually can happen can be shown by examples. 

#2 This was not so in our example I. 


SS a 


! 
i] 
‘ 
4 
‘ 
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positive s-values by termwise integration so that (6.6) certainly holds for suffi- 
ciently large positive s. Now with 2|A;| converging, (6.6) defines «(s) 
uniquely for all complex s (with singularities at the points s, and the points of 
accumulation of s,, if any). Since the analytic continuation is unique, it follows 
that (6.6) holds for all s. The series = | A, | must, of course, converge if (6.2) 
is to converge absolutely for ¢ = 0, and this proves the necessity of our condi- 
__ ¥(s) 
1 — ¢(s) 
verges, then w(s) is the Laplace transform of a function u(t) defined by (6.2), 
Since the Laplace transform is unique, u(t) is the solution of (1.1) by Theorem 2, 
The series (6.2) converges absolutely for ¢ > 0 since | Axe’ | < | Ax |e”! 
Finally (6.7) follows directly from (6.6). 

It is interesting to compare (6.7) with formulas (50) and (56) of Lotka’s 
paper [8]. Lotka considers the special case g(t) = f(t); in this case y(s,) = 


o(s:) = 1, and (6.7) reduces to A, = — 


tion. Conversely, if w(s) = is given by (6.6), and if =| A, | con- 


ie . If s, lies in the domain of con- 
¢' (sx) 


vergence of the integral g(s) = | e “‘f(t) dt, that is, if R(s,) > o then 
0 


(6.8) 2 de [ e* yf(t) dt, 

A; 0 
in accordance with Lotka’s result. However, (6.8) becomes meaningless for the 
roots with R(s,) < o, whereas (6.7) is applicable in all cases. 

Theorem 6 can easily be generalized to the case where the characteristic equa- 
tion has multiple roots. The expansion (6.6) (which reduces to the customary 
expansion into partial fractions whenever w(s) is meromorphic) is to be re 
placed by 


- At” Ay” pe ) 
(9) 0) = EAS G + goat teal 


where m, is the multiplicity of the root s, . This leads us formally to an 
expansion 





p7* 
(m, — 1)!)’ 
which now replaces (6.2). Generalizing Theorem 6 it is easy to formulate some 
simple conditions under which (6.11) will really represent a solution of (1.1). 
Other conditions which ensure that (6.9) is the transform of (6.10) are known 
from the general theory of Laplace transforms; such conditions usually use only 
function-theoretical properties of (4.9) and are applicable in particular when 


w(s) is meromorphic. We mention in particular a theorem of Churchill [17] 
which can be used for our purposes. 


(6.10) u(t) = 2 eae + A ‘ -... +e 


7. On the practical computation of the soiution. There are at hand two main 
methods for the practical computation of the solution of (1.1). One of them 
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has been developed by Lotka and consists of an approximate computation of a 
few coefficients in the series (6.2). The other method uses an expansion 


(7.1) u(t) = d walt 


where u,(t) represents the contribution of the nth “generation” and is defined 
by 


(7.2) uo(t) = g(t), Unq(t) = I Un(t — x)f(x) dx. 


Now the Laplace transform of wn+:(t) is y(s)e"(s), and hence (7.2) corresponds 
to the expansion 


_ vs) _ ~~ 
(7.3) Me 7= oe) ¥(s) Le (s). 
In practice the functions g(t) and f(t) are usually not known exactly. Fre- 
quently their values are obtained from some statistical material, so that only 


their integrals over some time units, e.g. years, are actually known or, in other 
words, only the values 


(n+1)8 (n+1)8 
(7.4) i ; [ fd, gn = : [ g(t) at 


are given, where 5 > 0 is a given constant. Ordinarily in such cases some 
theoretical forms (é.g. Pearson curves) are fitted to the empirical data and 
equation (1.1) is solved with these theoretical functions. Now such a procedure 
is sometimes not only very troublesome, but also somewhat arbitrary. Con- 
sider for example the limit of u(t) as t — ©; this asymptotic value is the main 
point of interest of the theory and all practical computations. However, as has 
been shown above, this limit depends only on the moments of the first two 
orders of f(t) and g(t), and, unless the fitting is done by the method of moments, 
the resulting value will depend on the special procedure of fitting. Accordingly 
it will sometimes happen that it is of advantage to use the empirical material 
as it is, and this can, at least in principle, always be done. 

If only the values (7.4) are used it is natural to consider f(t) and g(t) as step- 
functions defined by 


(7.5) KO) = fa, for ni <t < (n+ 138. 
g(t) = gn, 


In practice only a finite number among the f, and g, will be different from zero: 
accordingly the Laplace transforms y(s) and ¢(s) reduce to trigonometrical poly- 
nomials, so that the analytic study of w(s) = v(s) becomes particularly 


1 — ¢(s) 
simple. Lotka’s method can be applied directly in this case. 
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For a convenient computation of (7.1) it is better to return to the more genera] 
equation (1.3), instead of (1.1). The summatory functions F(t) and G(t) should 
not be defined by (1.2) in this case, but simply by 


(3] [3] 
(7.6) F(t) = do Sas G(t) = 2, n- 


It is readily seen that the solution U(t) of (1.3) can be written in the form 
U(t) = >~ U,(t), where 
n=0 


U,(t) = Gd), U,.(t) = [ U,(t — x) dF (zx); 
0 


in our case U,(t) will again be a step-function with jumps at the points ké, the 
corresponding saltus being 


k 
(k) (k) (k—r) 
uP=g, ui = Do Un ‘. 
Tr 


Thus we arrive at exactly the same result as would have been obtained if the 
integrals (7.2) had been computed, starting from (7.4), by the ordinary methods 
for numerical integration of tabulated functions. It is of interest to note that 
this method of approximate evaluation of the integrals (7.2) leads to the exact 
values of the renewal function of a population where all changes occur in a dis- 
continuous way at the end of time intervals of length 6 in such a way that each 
-change equals the mean value of the changes of the given population over the 
corresponding time interval. 
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ON THE JOINT DISTRIBUTION OF THE MEDIANS IN SAMPLES 
FROM A MULTIVARIATE POPULATION 


By A. M. Moop 
University of Texas 


It is well known [1] that in the case of a population having a single variate 
distributed according to a density function satisfying certain general conditions, 
the median of a sample is asymptotically normally distributed about the popula- 
tion median as a mean. It is the purpose of this paper to extend this result to 
populations involving more than one variate. Besides the theoretical interest 
of such a result, there may be some practical value in it when one is dealing with 
samples from a population for which the median is a more efficient statistic than 
the mean, as, for example, when the population variance is not finite. 

The complexity of the exact distribution of the sample median increases 
rapidly with the number of variates which describe the population; it is almost 
impossible to write out completely the distribution for the general case of k 
variates. For this reason the author has chosen to give first a detailed presenta- 
tion for the case of two variates, then use a condensed notation to establish the 
general result. This is a circuitous route, but it seems to be the only feasible one. 
A condensed notation is necessary for the general case, but presented alone it 
would be well-nigh incomprehensible. 


1. Distribution of the median in two dimensions. An extension of A. T. 
Craig’s [2] geometrical argument will be used to obtain the exact distribution of 
the sample median. Let us consider two variates x, and x2 with density function 
f(x: , X2) which shall satisfy the following conditions: 


1. f(t, %) >0 


+ [.5(% }) dx, = [_ fe,0) az, +0(t) 
mice n») dt, = [50 22) dz. +0(t) 


3. [ [se, Xe) dx, dxz = 1 


4. Each of the equations 


[i [ses 29 aes 


Nie 


1 
2 


[ [ f(xi, 22) dxidz, 


has a unique real root. 
268 
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If & and & are the respective roots of the two equations of this last condition 
then the point (£, , ) is defined to be the population median. It will be assumed 
in what follows that the coordinate system has been so chosen that é = 0 = &. 

Let a sample of 2n + 1 elements (ia , Tea)(a = 1, 2, --+ , 2m + 1) be drawn 
from this population. The sample median (%; , %2) will be defined as an element 
(not necessarily in the sample) whose 2, coordinate is the middle, with respect 
to magnitude, number of the set of numbers 2;. , and whose 22 coordinate is the 
middle number of the set of numbers 224. Now let us compute the probability 
that the sample median will lie in the rectangle 


Z:— $ dz < a < %; + 4 dz; 1=1,2. 


This rectangle will be denoted by R’”’. The remainder of the plane will be divided 
into eight other regions R, , --- , Ry as indicated by the dotted lines in Figure 1. 
The probability that an element will fall in the region R{” will be denoted by 


py - J Joc Seu 22) dx, dx, 


' 
' 
‘ 
' 
' 
' 
' 
‘ 
' 
' 
' 
‘ 
' 
' 
‘ 
' 
' 
Tr 
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- 


Fig. 


Neglecting terms involving differentials of higher order we have 


~~ = [ [, Hm, t2) dxedx, 


P2 = c [ S(ti, 2) dae dx, 


p' = | Fler, ¥) dey dh 


p” = f(%, %2) d%, dz. 
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We shall consider now that the sample is drawn from a multinomial population 
with probabilities p,, --- , p’’ and pick out those terms which give rise to 
sample median in R”. If the median is an element of the sample, then that 
element must fall in R’” and the other elements must fall in the regions R, , 
R., Rs, and R, in such a manner that 


M+m=n+m=n 
M+m=m+n=n 
or so that 
(2) m = ns and nm = 


‘where n; is the number of elements in R;. The probability that this occurs is 


2n 1 ! nh n n n 
(3) 2 Got? p' pi" paps ps 
ni rne=n . ° 





Now suppose the median is determined by two different elements of the sample, 
for example one in Rj and one in R;, then there must be n, elements in R,, 
nm, + 1 elements in R; , and nz elements in each of R; and R, with 

(4) mM+tm=n-— 1. 

The probability in this case is 

(5) Dips (2n + 1)! ny! ne, nitl, ne 


artetin-1 ml (m + 1)ime P2 P3 D4 e 


Continuing in this manner we obtain the distribution of the median, and letting 
D(H , Z2) represent the density function giving this distribution we have 


a a o ~ ” (2n + 1)! ny / ne 
D(&, %2) d%, diz = p"’ = me me (p~i Ds) (po Pa) 
i % o_# (2n + 1)! ny n 
(6) + (pspiP2 + PiPsPs) z m!(m + 1)!mP (pips) : (pep) ' 
(2n + 1)! 


+ (papi Pi + PsP2 Ds) 2 (pips) (p2p)”™. 


2. Asymptotic distribution of the median in two dimensions. As a simple 
notation 
A = Bil + O(1/Vn)) 
will be abbreviated to read 
(7) A =.- B, 
the dot after the equality sign indicating the omission of the factor 1 + O(1/+/n). 
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As is customary, the second term of this factor represents any function such that 
lim NO(1/N) = L < o. 
N-0 
In order to get an approximation to (6) for large n we shall use the normal 


approximation for the multinomial distribution and compute the sums (these 
cannot be put in finite form) by integration. We use then the well-known result 


@) [ot =. 14/020) exp (— 4S Aue) TH ae, 
II mm! 1 1 i 
where 
(9) zi = (m; — mpi)/~/m, 7=1,2,---,r-l1, 


(10) 


Returning to (6) it is to be noted that the fraction immediately following 2 
in the first sum has one more factor in the denominator than the corresponding 
fractions in the other sums. This first sum may therefore be neglected in the 
asymptotic form as it is of order 1/n in comparison with the others. We con- 
sider now the second sum in (6) and let it be represented by the letter S 

nitl ne 


ai ro? (2n _ 1) nine 
(11) S = 2n(2n + 1)pipe ll m! (m+1)1 mt”? P2° Ps Pas 


Employing (8) and omitting certain terms of order 1/n we have 
. 3 
(12) S =. 4n*pip; D |A/(2x)*} exp (- $ dX Ajj2i2; ) dz, dzz dzs, 


in which the A,;; are defined by (10) with r = 4, and 
(18) 2 = (nm — 2npi)/+/2n, i = 1,2,3. 
In view of the relations (2) between the n; we have 

a= V2n(t—-— mm) -—aA=u-A 

23= V2n (pi — ps) —%41=%— 4%, 


in which relations we have defined the new symbols u; and uz. It will be recalled 
that in (8) the factors dz; correspond to factors 1/+/m, we therefore let dz, and 
dz; in (12) cancel a factor 2n from the coefficient of the exponential, and after 
substituting (14) in (12) find that 


(14) 


_ tem{—afa(t+t4242) 
S =. Inpips E[A/(2r)'I exp | y[a(t+t+i+2 


2 2 2 
+2,(utm wy + tul ty ean, 
Ps P2 Pps Pa Pe Pps 


(15) 
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The summation can now be performed to within terms of order 1/+/n by inte- 
gration with respect to z, between the limits — © and + ©; this gives us 


, / ; 2 2 
g =. BP 4 /(L 4 bt 443) cxp{— ya two 44 
2r Pi D2 D3 Pp —4 Do 


+B— (Hee wy) /(Letyt ey] 
Ps Ps P2 Ps Pi P2 Ps Ps} 


At this point some new symbols are required. We let q; and q; represent the 
results of replacing Z; and 72 by zero in the integrals of the relations (1) 


(16) 


“a = [ , S(x1, X2) dx dxe qa = [ se. 0) da, 
qe = [. [ se, Le) dx; dre q = [ 10, Le) dx2 
(17) — : 
qs = [. [sa Le) dx, dre q3 = [ te, 0) dx; 
2 0 0 
“= [ [ fe, 2 dx; dx2 a= [_ 10, Le) dx, 
then 
(18) Q+R=~BtUH-UAtUN=H~HR+ GS =} 
and 
(19) H=%, g=%. 
Also we let 
(20) n= a+q, m=anta, 
(21) y= V2Qnaki, ys = V2n aks. 
We have now 
Di =: Qi; 2 = 1, 2, 3, 4, 
(22) Di =- gidks, i= 1,8, 
Dp =- gdh, i = 2,4. 


Also 
uy = /2n(4 — pi — pr») 


= Vin [se tq) dar, div, 























the 
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Vink | Slar, 08) da, 
23) ° 


. Vis, | fa, 0) da ; 


=. /2n ak, 
wii 





Similarly 
V/2n (pi — Ds) 
- —(y1 + y). 


The result of substituting (22), (23) and (24) in (16) with some further simplifica- 
tion using (18) and (19) is 


§ 
ll 


(24) 


2ngi qs ( 1 yi - 4(q — Q2)Y1Y2 + st) a 
25 8 am. —— shes OE Ct ae aninaeaninniinsiieneaimaninaes BUENa MTR 
(35) Qrv/ a9 ™ © 49192 — 


The other three sums of (6) will give rise to the same expression except that the 
factors 9192 will be different; it is clear then that 


ae | a e # 
D(:, Zs) d&, 0%, = - 2n(qige + gig2 + 9293 + 9394) 





2rvV 192 
1 yi — 4(q: — g2)yiy2 + #) " 
ii iy. i imamate stata ee dz 
x exp ( 5 = 1 dE. 
2na; 2 ( ai Fi — 4(q: — g2)Q1a2%1%2 + #3) n 
26) =. SO) ee i PY 
; Van 49192 ; 


(27) 


1 ( Lyi — 4(q — 92)yiy2 + =) 
. —L exp ( -1 4 — 4 — Gh +H) 4, 4... 
Qar-V/ng - 491 q2 —_ 


This is the asymptotic form for the distribution of the median in two dimensions. 


3. Distribution of the median in k dimensions. We consider now a population 
characterized by a density function f(7 , --- , zx) defined over a euclidean space 
of k dimensions satisfying conditions like those required of f(z; , x2) in section 1, 
and we assume that the population median is at the origin so that the integral 
of the density function over any half-space determined by a coordinate hyper- 
plane is 4. 

A sample of 2n + 1 elements will have a median (%,, --- , %) each coordinate 
of which is the middle number of the set of numbers giving the corresponding 
coordinate of the elements of the sample. To obtain the probability that the 
sample median lies in the hyperparallopiped %. — $d%a < ta < Za + } dz. 
(« = 1, 2,--- , k), we divide the space into 3° regions by means of hyperplanes 
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perpendicular to the coordinate axes through the points Z, + 4} d%, on the co- 
ordinate axes. These regions are illustrated in Figure 2 for the case of three 
dimensions. The coordinate axes have been omitted in this figure. There 
will be 2° primary regions denoted by R,, R2, --- , Ree corresponding to the 
octants of the figure; k2‘~* regions with one differential dimension denoted by 


Ri, Rz, --+, Risk-1 corresponding to the quarter slabs of the figure; (5) ii 


regions with two differential dimensions corresponding to the half strips of the 
figure, and so forth. Probabilities associated with these regions are defined by 


i) 
pi? = Jf, Har, +++, m0) dan +++ da. 
Ry 





If the sample median is determined by k different elements of the sample there 
will be one of these & elements in each of k regions R; whose differential dimen- 
sions are mutually orthogonal and the other elements of the sample will fall in 
the regions R; in such a way that n elements of the sample will lie on either 
side of any of the k hyperplanes z, = %,. The probability of this occurrence 
for a particular choice of k of the regions R; is 


k / 2k 
(28) s= Tle. Det TT pe 
a= Il Nn! i=l 


in which the 2" indices n; are subject to k independent restrictions of the type 


(29) wn =n — Ca; 
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where C. is an integer such that 0 < ca < k, and the prime on & indicates that 
the sum is to be taken over all n; on one side of a hyperplane 7. = %,. 1; is 
the number of elements in FR; and besides the k conditions (29) we have also 


ak : 
(30) din = In-k +1. 


In order to include all ways in which the median is determined by k different 
elements of the sample we must add together 2*“~” sums of the type (28). If 
the median is determined by less than k elements, say k — h elements, then the 
fraction (2n + 1)!/IIn,! will have h extra factors in the denominator and hence 
the sum will be of order 1/n” as compared with that of (28) and may be neglected 
in obtaining an asymptotic expression. 

Thus we need only find the limiting form of (28) 


S = (Qn + 1)(2n) --- (Qn — k + 2) IL». >> ee + DTT pt 


which after substituting (8) and neglecting terms of lower order becomes 
2k1 
(31) S =-(2n)* TT pi. > (A/(2e)™)! exp (—4 DY Azzizs) I] dz;, 
in which the A;; are defined by (10) with r = 2* and 
(32) 2; = (nx — 2npi)/~V/2n, i=1,2,.--,2—1. 
Now we define 
(33) Ua = V 2n($ - 2’pi), a=1, 2, es k, 


the 2’ having the same significance as in (29). These conditions (29) may now 
be put in the form 


Za =-Ua — L,(2), 


in which L,(z) is a sum of a certain subset of the variables zi41, --- , Z#-1. 
Care must be taken in labeling the regions R; in order to be able to solve for 
a4 ,---, 2 in this form. After substituting these relations in (31) we replace 
k 


II dz, by (1/2n)*” and perform the summation to within terms of order 1/+/n 
by integrating the remaining z; from — © to +; the result is 


k k 
(34) S = -(2n/2n)** TT pi, VB exp (-3 . Busta), 


in which the Bag are functions of the p;, and B = | Bag|. As in (17) and (20) 
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we define 
, «=| flm,i-,a)Mdr, 


(35) qi = | f(t, +++, %)M' dre 


a, = / f(ai, +++, &)II' de. = 2'g;, 
Zqm0 


in which R; is the set of regions bounded by the coordinate hyperplanes R; 
are regions into which the coordinate hyperplanes are divided by the remaining 
coordinate hyperplanes. II’ indicates that one of the differentials is omitted 
and the variate corresponding to that differential is put equal to zero in 
f(a, +++, 2%); 2’ indicates the sum over all g’ determined by regions lying in 
the hyperplane zz = 0. It is clear that 





Pi ™ Gi 
an IT pi. =- IT 9.4% 

Uq =-v/2n . Sag Agha = >, dap Ye, 
where 


bag = +1 or 0, and Ys = a/ 2nasks . 
Making these substitutions in (34) we have 
k k 
(37) S =.(2n/2n)*” TT ai, VC exp (=n Z Costes) I] azz, 
1 1 


and adding together all possible sums of the type (28) we have the asymptotic 
form of the distribution of the sample median 


(38) D(&,-+-+, %) I] dz 
. (2n/2m)*” II dar/C exp (=n > Casas) Il dia 
(39) -(1/2m)*?4/E exp (—4 d) Casyays) [] dye, 


in which the C.s are functions of the q; . 


4. The case of three dimensions. - The computation of the coefficients Cas 
of (39) requires the evaluation of a determinant of order 2° — k for each one of 
them. This work was quite laborious even for k = 3 and the author made 
no attempt to find their explicit expression for larger values of k. 

If we let a subscript + indicate integration of the density function 
f(x: , x2, 23) from 0 to ©, and a subscript—indicate integration from — © to 0, 





es R; 
Lining 
nitted 
ro in 
ng in 


ptotic 


I dz. 


ts Cag 
one of 


made 


netion 
to 0, 
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as for example, 


SuH- = [ [ [teu 1, 2) dndrdn, 


then the q; of (35) will be defined as follows 


M = fri qs = fie 

Qe = fxy- ge = fi 
(40) 

@ = fr @ = f—+ 

Mw = fx— % = f__ 


The coefficients C.g may be written 
DC = 2(q1 + 95)(Ge + Qe) 


DC = 2(q1 + Qs)(G2 + Qa) 
DCs = 2(q: + 2)(gs + 4) 
41) DCw = 9395 + U9 — 197 — 429s 
DCs = G95 + 9497 — 19s — 9398 
DC23 = 293 + aay — MQ, — 9598, 
where 
D = crasanau( 2 + - + - + 4) + asavaran(2 + : + i + 1) 
+ 2(qs + 96)(q7 + 98)(Q192 + 9594) 
(42) 


+ 2(qs + 97)(G6 + 9s)(qigs + 929) 
+ 2(gs + as)(@6 + 97)(qige + 929s) 
+ 8(q1949697 + 92989598) 


(41) and (42) can of course be put in different forms by using the four relations 
between the g;. The a. of (38) are defined in (35); fork = 3 they are 


aq = [ [. S(O, x2, 2X3) dx2dxz 
(43) Oy = f. [ fles, 0, 22) drida, 


a3 = [ [ se, Xe, 0) dx de. 
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5. The normal distribution in two dimensions. If the density function of 
the second section of the paper is normal 


(44) f(a, 22) _ 1/(2ro1020/1 — p*) exp| 1 (= ie 20 Vie + )), 


— 20 =p) ae * Gh 
we find that the parameters of (26) are 
ae Pe _1 er 
m= 74 + >5- sn P; | = P; 
(45) : 1 


ay 


=> —— - a= - 
/ 29 0’ V 21 oe 


These give an interesting result—the correlation coefficient of the asymptotic 
distribution of the sample medians is 








(46) tas = des 
T 

hence 

(47) |em| S |p| 


the equality sign holding only when p = 0 or +1. 
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SAMPLES FROM TWO BIVARIATE NORMAL POPULATIONS’ 
By Cuune Ts1 Hsu 
Columbia University 


1. Introduction. In multivariate analysis involving p variates, or in analysis 
of variance of m samples from univariate populations, we are often interested 
in the hypothesis of the equality of variances; viz., that 


O1 = 02 = +++ = Gy, in the case of p variates; 


1 = 02 = +--+ = Om, in the case of m samples. 


As a matter of fact, it seldom occurs that these hypotheses are true, but the 
ratio between the variances might be known. 
Hotelling [5] has suggested that if 
oi/ky = 02/ky = --- = on/km = 0", 
where the k’s are known constants, we can apply the transformation 


, 
= Wit, 


where 
wV/k; = 
so that after transformation the variances become equal, i.e., 


eS , 
O1 = 62 = +++ = Gm, 


and the required analysis can be carried out. This method is similarly ap- 
plicable in the multivariate case. 


In a previous paper [7], I developed a series of hypotheses concerning samples 
from a bivariate normal population under the assumption that 


O1 = 02. 


In case o1/k; = o2/k2, where k, and k, are two distinct known constants, 
similar results may be obtained by the use of the transformation t= Wit; 


Zy = Wet, ; where wiv/k, = We/ky = 


1 Presented to the american Mathematical Society at Washington, D. C., May 3, 1941. 
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In multivariate analysis, the hypotheses usually of interest concerning correla- 
tion coefficients may be classified in two categories, viz., 

(i) that the correlation coefficient is equal to a specified value, e.g., in 
simple correlation pi2 = po , in partial correlation, pi2.3 = po , in multiple 
correlation, pi.3 = po, or in correlation between two sets of variates 
[4]’, Q = Qo; of special interest is the hypothesis of the vanishing of 
such correlation coefficients. 

(ii) that two given correlation coefficients are equal, e.g., (1) correlation 
coefficients p; and ps in the correlation matrix of a multivariate distribu- 
tion are equal (Hotelling [6]), or (2) the correlation coefficients py. and 
pi2 in two bivariate populations are equal. 

R. A. Fisher in his earlier paper [3] introduced the transformation z = 

” 
, log —! which provides a very satisfactory, though approximate, method for 
the comparison of two correlation coefficients. Brander [1] treated the same 
problem by the method of the likelihood ratio criterion. 

The present paper is an attempt to obtain different criteria by the likelihood 
ratio method (Neyman and Pearson [9], [10], [11]) for testing, by means of 
samples, the equality of correlation coefficients in two bivariate normal popula- 
tions under the following sets of conditions: (1) 0: = o2 and 0; = 02 ; (2) 0, = 7, 
t, = £ and oj = o2, ; = &. The results may be extended to the cases (3) 
01/k: = 03/kz and o3°/ki = o9°/ke ; (4) o1/k: = o3/ke, &1/k: = &/k2 and o'/ky = 
o2'/ke , &1°/ki = &’/ko, where the k’s are known constants. 





2. The hypotheses. Two samples, each being of two variates (2 , 22) and 
(x; , x2), of size N and N’, are supposed to be drawn at random, respectively, 
from two independent normal bivariate populations, with the following distri- 
butions: 


eeeenllimeend = \- 1 (2 - =) 
Qroioe/1 — p 2(1 — p*) 1 
= (2=2)(258) (058) 
01 02 02 
satrap tatlls) 
(2) 2ro1020/1 — p? 2(1 — p”) o1 


=, S — &) (* = " + CG [ =I}, 
01 02 02 


where £1, &, 01, 02, p; £1, &, 01, 02, p’ are the unknown parameters of the 
populations. 

The hypotheses to be considered in the present paper are: 

H,: Assuming o; = o2 and o; = az, to test p = p’. 

H.: Assuming o; = o2, & = &, and o; = 02, & = &, to test p = p’. 





(1) 


2 See bibliography at the end of the paper. 
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) and 
ively, 
distri- 


NORMAL POPULATIONS 281 


The derivation and the distribution of the criteria for testing these hypotheses 
may be simplified by the following simultaneous transformations: 


(3) X = zi (x1 — 22) “v; (x + Ze) 
(4) X' = - (xi — 22) Y' = = (xi + 22) 


The corresponding normal bivariate distributions in the transformed variables 
(X, Y) and (X’, Y’) are obtained, viz. 





en ma lS) 

QnoxoyV 1 — pxy 2(1 — px) ox 
_ -_ — 

too) 54) «(%S) Dace 

ox Oy Oy J 

STE wills) 

, » 7 exp a 12 4 
QroxoyV 1 — eae 2(1 — pxy Ox 


, nie , ‘aie , , —_ *\2 
aio 5859) + ("GAY aver 
Ox Cy Oy 


The conditions corresponding to 


(5) 





(6) 











(7) o1 =o: and o, = 02, 
are that 
(8) pry =0 and pry = 0. 
Also, for a given p and p’, we have from (7) 
(9) oy = yor and oy = ‘ox, 
where 
(10) ya pth and y= Te, 


Following the notation of (9) and (10), the hypotheses H; and H; corresponding 
to H 1 and H, are: 


Ay : Assuming pry = 0, and pry = 0, to test 7 = 
H;: Assuming pry = 0, — = 0, and pxy = 0,’ = 0 a 


3. The derivation of the criteria. Let (x1; , 22:)(21; , 22;) be the measurements 
of the characters on the 7th and jth individuals in the two samples from their 
respective populations. After transformation, the corresponding measurements 
become (X;, Y;) and (X;, Y;). Let p(E) denote the joint elementary proba- 


oseas ae —— 


Se 


wre tee ete Se ee ee 
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bility law of the N and N’ observations, H = (Xi,---,Xw, Yi1,-++,Yy; 
Bag ++ ple, OH, +++ Fp. 

Following Neyman and Pearson, we shall use 2 to designate the class of ad- 
missible populations under conditions which can be assumed to be satisfied in 
any case; and w to designate a subclass of 2 under conditions which are satisfied 
only if the hypothesis to be tested is true. 

Thus for H’, & specifies for pry = pxy = 0, any real values of &, n, &’,-n' and 
any positive values of ox, oy, ox, cy ; w specifies pry = pxy = O, any real 
values of &, 7, &’, n’ and any positive values of cy and y which are defined by (9). 
While for H’, Q specifies pry = pry = 0, & = & = O, any real values of 7 and 7! 
and any positive values of ox , oy, ox , vy ; specifies pry = pry = 0, =#’ = 0, 
any real values of 7 and n’, and any positive values of cy and y which are defined 
by (9). 

For our hypothesis H; , the values of the parameters required to make p(Q) 
@ maximum are: 


=X, f= 
fs, « 


iy 


A * 
’ Gx = &z, Cy = 8y 


/ 


, 8y. 


wet 


af _ “ip 
’ Ox = &, Cy 


Thus p(2 max) = (ty ae 
2x Gee & 
To obtain p(w max), let us define, according to the notation in the writer’s 
previous paper [7], 
2Ys; 82 ’ 2Y’s: 83 


‘+s ” ‘ap bap 





and 


1—R, > 1-R 


8x 





i= 


@ a 
bq to] ne 89 


Then the values making p(w) a maximum are: 
=X, 4 =Y, of = 38x4+ 4) 
f= 2%, =’, oy = ter + v’) 
and ¥ is the positive root of the equation 


(N + Ny’ — (N - N’)(u — u')\y — (N + N’)uu’ = 0 


or 
g = NON Yu — wv) + VW — NY — WY + AN + Nun! 
(11) 2(N + N’) 


= V1; say. 


’ Yy ; 
of ad- 
fied in 
itisfied 


-n’ and 
ny real 
by (9). 
and 
t’ =0, 
defined 


ce p(Q) 
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Then 
N+N’ N — Nn’ 
p(w max) = (2) [_ 2, | | avin 2 | gunn 
2a (v1 + u)sx (m1 + u’)sx 
and the likelihood ratio criterion for the hypothesis H; is 


_ p(w max) _ | aye, 211 Sr | 20/71 8y avn 7 


p(Q max) (yi + u)sx (yi + u (yi + u’)sx 


2 EF nel E =—y 
itu m+ u’ 


For H; , the values the parameters to make p(w) a maximum are: 


(12) 


2 
4 = ; g2 = 1 ox Gy = 8y 


N 


1 \*+¥’ / NN’ a 
pome) = (5) SEE 


Similarly, if we write 


_ 2y8ise — $(#1 — 4)” Ri = 2y'si8a — HG — a) 
si + 83 + 4(t1 — %)”’ si + 83° + 4H — 42)” 


1+R, ,_Nsy _1+R 
~ oie ick’ 


the values to make p(w) a maximum are: 


q= Y, oy = ne + v) 


=P, oF = 5 2X" +0) 


_ (N — Nv — 0’) + V(N — Nv — v’)? + 4(N + Nu’ 
2(N + N’) 


= 2, Say. 


_(1V"T INV 9 i 2N' V2 : 
aa (=) |Z + 5X2] LO + vex? ’ 
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and the likelihood ratio criterion for the hypothesis H2 is 
A» = p(w max) | 2\/N728y ial 2v/ N’ ¥28y : 








(14) ~ p(@max) ~ LG + o)V/ axe] Ls + vVER" 
ene 
oe +v ¥2 + v’ . 


The caseN = N’. The above criteria \; and d- cannot in general be expressed 
simply, but when N = N’, by (11) and (13) 


n=Vu', n= Vo’, 





and 
it | AY] ae | al 
Vat Vurd ’ Wot vord’ 
or we may express as monotonic functions of \; and dz, 
(15) Ly = vOtw) = UY = 4 es 
eV) 
u u 
(16) Ig = MM" = 7 * 


V+ 5) 


Thus, 2’s, L’s, or their functions = “3 may be used as the criteria in the 


present case. 
Furthermore, if we introduce, 


(17) z2=tlogu, and 2’ = }logw’, 


we have 


/ / 


U 4 ’ 
4(z— 2’) =hlog— or zee". 
u u 


Thus LZ; can be written in terms of z and z’ 
(18) Ly = 4/(e**? 4 et”) = 1/cosh? #(z — z’) = sech’ #(z — 2’), 
and z — z’ = w, say, may be used also as a criterion for H,. 


We shall now proceed to obtain the distributions of some of these statistics. 


4. The distributions of u/u’ and v/v’. Since Ns?/o? and Ns*%/o% have inde- 
pendently the x’ distribution with N — 1 degrees of freedom, 
Sy oy x2 pe 1x2 


2 2 2 
8x OxX1 X1 


and u/y has the F distribution with degrees of freedom f; = N — 1, fe = N — 1. 
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Similarly, u’/y’ = xs’/x:° has the F distribution with the same numbers of 
degrees of freedom (since N = N’, in the present case). 
If the hypothesis H; is true (i.e., y = 7’) 


(19) u _xaxt _ O16 _ 2 


ul xix.” 010, 


where 0;(—3x?) or 6; is distributed as 
1 gett ap. 
(20) iva) 6s ee * db, 
with a; = 4(N — 1), and 2,(= 6162), ze(= 6162) follow independently the Wilks’ 
z-distribution, [14], which we shall study in detail for the present case. 
Distribution of z when p = 2: Consider 


z2 = Bo,O.--- 05. 


Wilks has succeeded in integrating the distribution of z for the case p = 2 for 
special values of a’s, e.g., a1 = $(N — 1), ag = 3(N — 2). Now we want the 
distribution of z when p = 2 and for any values of a, and then for a; = a@ = 
4(N — 1). 

By (20) the joint distribution of 6, and 6: is 


1 G1—1 —6; pa2—1 —8 
——__——~ 6" e *' 63? * 0; d0.. 
T(ai)I' (an) 1 2 @€ 1002 


Applying the transformation z = Bé,6. , 1. = 6,, the joint distribution of 1, , z is 


1 @;—1 —v ( 2 , —2/Bv} dv; dz 
——_——. pe — e =. 
T'(a;)T' (a2) Buy By 
Integrating v, from v,; = 0 tov, = ©, we have the distribution of z, viz., 


23" dz - @4—ag—1 —v1—2/Bv, 
(21) BT (ai)T(as) : v é dv. 


In order to evaluate the integral of (20), consider the transformation », = y’, 
dy, = 2y dy, we have 


(22) Io = 2 | pent oP dy. 
To evaluate I for any a’s, by putting y = 1/2, dy = —dz/z’, we have 
© —sx2/B—1/z2 


é 
(23) Io _ 2 | pea dz. 


Consider 


(24) T'(a, — a + 4) ins f ere yaad dy. 


x (a@,—a2) +1 0 
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Then 


Io TV (a1 — ae + 3) = 2 I ge ee? /B+1/z9) dx I ev yaa dy 


a [ yt dy [ eg lGlBtyst+1/23) dz 
0 0 


a a dy 
= vz | ge rete pment —— 
0 V2/B+y 
Since by the substitution i +y= Vi +yory=2+2 5% dy = 


2(z o i 5) dz and therefore 


IoT (a1 — ae as 3) = avs | ew #/ B+2) (= + or oe dz, 
0 


aww * uv aTBt2) (« 4 2r (i) ae 
(25) fo= ae B 


Hence, z is distributed as 
Q4/a gt evils « ais ( of; : il 
26 —— “13 + 2 panne" de. 
= B” T(a,)0'(a2)T (a1 — ae + 4) I , B wl . , 


We infer from this distribution that when 2(a,; — a), ie., the difference of 
degrees of freedom, is odd, the integral can be expressed as a terminated series; 
but for even values of 2(a: — ae), the series is infinite. 





When B = oa = }(N — 1), a = }(N — 2), (26) is reduced to 
(27) /7 A” 22 te? Az 
I'(a;)T'(a2) : 


which is Wilks’ é distribution, [15], for p = 2. 
When B = 1 and a, = a@ = 3(N — 1), it becomes 


(28) Jae tve | ” MAE + 2) tae, 
T'(ai)T' (ae) 0 


which is the distribution of z involved in (19). 

Since (28) can apparently not be simplified, I have been unable thus far to 
find in manageable form the distribution of the ratio z;/z, and therefore of u/u’ 
in this case. However, it would be simpler to use the alternative criterion 
w = z — 2’ for the hypothesis H,. The distribution of w will be taken up in a 
later section. 
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The distribution of v/v’: Since Ns;/oy and 2X°*/cz have independently the x’ 
distribution with N — 1 and N degrees of freedom respectively, therefore, 


xxi 
and — s/w —— ! has the F-distribution with f, = N — 1 degrees of freedom and 
h= ’ 
Similarly - / ~_— has the F-distribution with degrees of freedom fi and fe 


as above. 
If the hypothesis H; is true (i.e., y = vy’), 


where each @; is distributed as in (19), but with a, = 4N and a = 3(N — 1). 
We can infer from (27) that th = 4./z and } = 4+/z have independently 
the x’-distribution each with 4a, or 2(N — 1) degrees of freedom, and t/t = 
/ 1/2 = +/v/v' follows the F-distribution with degrees of freedom f; = fo = 
2(N — 1). The 5% and 1% points of the F = v/v’ may be obtained from 
_Snedecor’s table ({12], p. 174). 


5. The distribution of y = log z. Wald [13] has suggested that the distribu- 
tion of z = B6,6. ----6, for any a,’s (¢ = 1, --- , p) may also be obtained in- 
directly with the aid of the characteristic function. A similar method has been 
applied in a recent paper by Wald and Brookner [14]. Consider the trans- 
formation 


(29) y = logt = log Bé,6.--- 6,. 
The characteristic function of y is 
g(t) = E(e”) = E{ (B06 --- 0,)°} 
(30) _ BYT(a + Ola + - - Tap + 1) 
T'(a;)T'(ae) - Fa) 
Thus the distribution f(y) dy is given »? 
a) sw=s [eran df wer YE atoy 


Without loss of, generality, we may ‘tie a = &® 2--- 2 a, > O and let 
a+ t= —?’, then 


eptio ’ ’ P 
(32) fly) = 2. eB T] ra; — a, — t’) dt’, 
21 i=1 


Ap—ico 


Pp 
where c, = e?”B™* / |] Tr(a,). 
t=l1 
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The integration can be carried out by the method of residue along the contour 
C, bounded by the line x = —a, and that part of the circle with center at 
origin and radius 7, which lies to the right of the line x = —a,. The integral 
of the function e’”B~“ [][?_, r'(a; — a, — t’) along the arc converges to zero 
as the radius of the circle tends to infinity (Kullback, [8]). Hence the integrals 
along the vertical line x + a, = 0 and along the closed contour C are equal. 
Then we may write 


— yt’ p-t’ . el: = ’ 
(33) fy) = - <2. fe BT] r(a: — a, — tar, 


and its value is c, times the sum of the residues at the poles within the con- 
tour C. 
For the present purpose, p = 2, we have 





—agtt 


(34) fy = 2. era — @ — e)r(-2) ae’. 


212 — 49+400 


We shall study the integral of (34) in more detail in the following cases: 
(i) a, — a = 3. By the duplication formula 


ry — ¢)r(—t’) = 2°" 4/xT(—2¢), 


and the function 


— NIN 
r(—2t’) = lim (oars 1) (a EN)’ 


has simple poles at the points 0, 3, 1, 3/2,---. The residue at t’ = m/2, 
where m is zero or a positive integer, is (—1)"""/2-m! and (34) becomes 


fly) = Vie(1 — 2d 4 Later — Lotgn 4 ++) 


=~ /r oe 


The distribution of z = e” is 


(35) 


, Qn/ motte tv 
2 eieieiinatieonn 
(27 bis) T(ai)F(an) 
(ii) a, — & = m+ 3. The function 


(a, — a — t’)T(-t’) = (m— 4 —-t)(m—§ -7¢/)--- 9 — “TG — “*)T(-#) 
= 2 /x(m — 4 —t)(m—§—t/) --- (4 — Y)T(—2¢) 
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has simple poles at 0, m, m + 4,m + 1,---,and 


_ (2m — 1)! y\m 1 a P 
fy) = V5 Ce laa DR (2’e”)” + + oom +1) (2’e”) +4 


1 y\m+ 
7 Om + aya?” _* | 


= Via| oma = - (2e et, 


2™—(m—1)! 2" < (Qm : y)7! 


This agrees with the expansion of (26) when we put a; — a, — } = m. 
(iii) a3 — a@ = 0. The function 


ne (N1?N 
TOON (88 (ONE AF PEN 


has poles of the second order at the points 0, 1, 2, 3, --- and 
© d , 
fy) = % De Ae — ye M(t) P my 
y=0 
(iv) a, — a@& = m. The function 


Tim — t’)T(-t’) = (m — 1 — t’)\(m — 2 — #’+).-- (1 — #)(-#)[T(-4)F 


has finite simple poles at 1, 2, --- ,m — 1 and poles of the second order at m, 
m+ 1,---,and 


m—1 


fly) = ¢ Dy {(t — ve" T(m = t)P(—t)}rny 
+ ce. z {SK ‘— ye *T(m — eon(—e)} 


6. The distribution of w = z — z’ or ¥ = cosh w. Since the distribution of u 
is given in [7] as 


1 u 4N-3 u —(N-) 
(39) WU — 1), KV — DI () (1 a *) ™ 


therefore, by ener (17), we have that the distribution of z for a given 
f= 4h logy = } log iis 


(40) 


as sech” (z — ¢) dz, 
n 
B(3,%) 











290 CHUNG TSI HSU 


where n = N — 1. The distribution of z has been given by R. A. Fisher [3] 
for n = 1 and by Delury [2]. Similarly, the distribution of z’ for a given ¢’ ig 


1 ' 
(41) ———~ sech”™ (z’ — ¢’) dz’, 
B(3 *) 
2’ 2 
where n’ = N’ — 1. 
In case-n = n’, the joint distribution of z and z’ for a given common ¢ is 


C dz dz’ 


(42) Csoch"  — $) noch! G' ~ dads’ = SoG eal @ — 8) 


where 1/C = E (G. a 


By the transformation Z = }(z + 2’), w = z — z’, we have the joint distri- 
bution of Z and w, 


C dzdw 2" C dz dw 


*) [cosh* (z — ¢) cosh (z’ — ¢)] - [cosh 2(2 — ¢) + cosh w]"” 
Integrating with respect to Z from — © to ©, we have 


n dz 
e dw | [cosh 2(2 — ¢) + cosh w]* 


(44) _ . 2dz 
“5 Cee lo [cosh 2(2 — ¢) + cosh w}* 


= 2"C dwl,, say. 


Applying the transformation ¢ = 2(2 — £), ¥ = cosh w, the integral of (34) 
becomes 


aa de 
1, = | (cosh ¢ + y)”" 


Substituting cosh ¢ + y = ey, we have 


I, = 


- oy cece rer 


ry — 9) rw _ STi 0) dé. 


Comparing (35) with the hypergeometric function 


I'(b)T'(c — b) 
I'(c) 


-ginl " 


(464) J= F e411 — 0)? (1 — 2) * dé = F(a, b, c, x), 
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we have b = n,c — b = 3, a = 3, and therefore (35) can be expressed in terms 
of a hypergeometric series as 

_ T(n)r) 1 ( ), 
n= Ta +h) WIP” inn th TS 


= 
oa 


(47) 


The series (37) is convergent since is less than unity. Thus the distri- 


bution of w, from (34), is 


2” CT (n)T(4) 1 cosh w — 1 
(48) hath Galetir’ (, a9 + RST a 


and the distribution of ¥ = cosh w is 


2"** CT (n)T(4) 1 y- ) 
Tin +) eee cin (hn + 4 Sot) ay. 


We notice that the distribution of y expressed in (39) is very similar to the 
r-distribution expressed in terms of hypergeometric series, except that in the 
first case the argument is ona , while in the second case it is iz where 
p = pr. Hotelling [5] has obtained a very rapidly convergent hypergeometric 
series for the distribution of the correlation coefficient since |p| < 1. But 
for the distribution of y, we cannot obtain a more rapidly convergent series than 
(39), since the values-of y lie between 1 and o. 


(49) 


7. Summary and remark. Two hypotheses concerning the comparison of 
correlation coefficients of two samples from bivariate normal populations have 
been considered. The appropriate test criteria for each hypothesis have been 
derived by the use of a transformation of the variates. The distributions of 
certain of the criteria have been obtained in the special case where N = N’. 
Incidentally the distribution of Wilks’ z for p = 2 and any values of a; and a, 
has been derived. 

Again though we assume throughout the paper that 0, = o2 and o; = o2, the 
tests can be generalized to fit the case where the ratios 01/02 = k, 01/02 = k’ 


are known, but are different from unity. In the latter case we can apply the 
transformation 


Yi = WN, Yo = Wete ; 
vi = with ’ y= Ws ; , 
where 
wiki = wke = 1, wiki = weky = 1, 


80 that after transformation the variances of each pair of y’s are equal. 
The writer is deeply indebted to Professor Harold Hotelling and Dr. Abraham 
Wald for their advice and suggestions in the preparation of this paper. 
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ON RANDOMNESS IN ORDERED SEQUENCES 


By L. C. Youne 
Westinghouse Electric and Manufacturing Company 


It is frequently desirable to examine an ordered sequence of measurements 
for the presence of non-random variability, concern over any particular type of 
variability being limited. Unless the sequence is one containing replicated 
observations, current methods of analysis often restrict an investigation to 
tests for specific forms of variability, such as particular orders of regression and 
periodicity. In order to simulate replication, arbitrary grouping of data is 
occasionally used and followed by some test of variance; this practice, however, 
is likely to add an element of bias to the investigation. 

Under these conditions, it would be convenient to have the means of testing a 
series for the presence of general regression, before proceeding to test for that of 
a specific type. It is the purpose of this paper to present, as briefly as possible, 
a statistic designed for this preliminary type of examination, and to demonstrate 
its application. 

If a given sequence of measurements be denoted by 


X,,X2,°-++, Xn 
then the magnitude of 


n—1l 
» (X; — Xin)” 


i- n ’ 
2 >> (X: — X) 
1 


will be dependent upon the arrangement of the n observations upon which it is 
based. C will have n! possible values for a given sample, corresponding to the 
number of permutations of n items. 


1. Moments of the distribution of C in terms of the moments of a 
finite sequence. Writing C in terms 2,--- ,2,, representing the devia- 
tions of X,, --- , X, from their sample mean of m measurements, 


a—l 


7 (x; = i41)” 


1 oe 1 
2 >> 2? 
1 
n—l 
ait+a2,+2 De times: 


22a 


1 
293 
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In order to find the mean value of C for a given sample, it must be summed 
over all values obtained from the n! permutations of the measurements. 
Dealing with the numerator alone of the expression given above: 


n—l n—1 
as E +2,+2 » ru | = Dietit Leta t+2 Qu» » Us Tiss, 


where 2, denotes summation over the n! permutations. 

There are n values of z;, and n! arrangements. Each value 2; is 2 in 
(n — 1)! of the arrangements: the same reasoning applies toz,. The first two 
terms of the summation, therefore, will be 


Dti = Dt = (n— 1)! dizi. 


With regard to the third term, there are 2(n — 1) of such cross-products for 
each arrangement. Since the summation is taken over n! arrangements, 2 jr; 
will be different than x,r;, and should be considered a separate term. Each 
2(n!)(n — 1) 

n(n — 1) 
arrangements, since there are n(n — 1) possible cross-products among n different 
items. The third term, then, will be 


n—1 n n—l n 
Ss. (= nz) = 2(n—1)!>> Doan = —2(n — 1)! Dd 2, 
1 1 1 1 
from which it may be seen that the mean value of C is zero for any sample. 


The same method may be applied in order to find the second and higher 
moments of C. Squaring the numerator of the expression and expanding, 


n—l 2 
r,| 2! +2,4+2> nts | 
1 
n—l n—l n—1 2 
= x, | z! +28 + 2zriz%, + 42 Do rizins + 42%, » Li tins + (5 nits) | 
1 1 
Performing the summation 2, term by term we obtain 


n—l 2 n 2 n 
ke i +2,4+2> nites 2(2n — a> 2) —2n > xi 
eidllidcneenncesiadaiadenaTlamcssaasilll di semenitensaiitlatseiill aititioeaninttenae® 


crossproduct term, therefore, must occur times throughout the n! 


1 
n! n(n — 1) 
whence the second moment of C for any sample is given by 
2n — 3 — m/m3 
2n(n—1) ’ 
where m, and m, are the second and fourth moments, respectively, of the n 
observations about their mean. 


In like manner, the third and fourth moments of the distribution of C for a 
given sample of n observations are found to be 


M, = 
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~ 6 +4(n — 3) M4 9™ _ 3 ™ 
M,=_____— 2 mms 
4n(n — 1)(n — 2) 
1 


2 
We 0 nti, Gin OP ~ eln ~ 
: NCES HESS | eo a me . 


’ 


— 24n(3n? — 17n + 27) ~ = a + (8n*° — 45n? — 23n + 210) ~~ 


+ 16(2n® + 5n — 21) “8° + 4(17n® — 37n + 42) 
Me 


— (7n’ + 13n — 6) ns). 
me 


2. Distribution of C for samples drawn from a normal universe. The 
first four moments of the distribution of C for samples drawn from a given popu- 
lation may be derived from the above formulae by substituting the mean values 


of ms =k etc. of samples from such a population. For normal samples con- 
ms Me 
taining n observations, for example, the following mean values apply, as obtained 


by the method presented by R. A. Fisher [1, 2]: 
m5 6(n _ 2) 


mi (n+ 1)(n+3)’ 

ms _ 3(n — 1) 

ms (n+ 1)’ 

my _ 3(3n + 23n* — 63n + 45) 

ms (n+ 1)(n+ 3)(n+ 5) ’ 

msms _ 60(n — 1)(n — 2) 
mi (n+ 1)(n + 3)(n + 5)’ 

ms _ _15(n — 1)” 

m; (n+1)(n+3)’ 

ms _ 105(n — 1)° 

m, (n+ 1)(n+3)(n+ 5)" 

Replacement of the sample moment ratios by the mean values of those ratios 

for normal samples yields the following moments of C: 


n—2 
M, = 0, M, > (n—1)n +1)’ M; = 0, 


3(n® + 2n — 12) 


Ms Get Det DetD: 
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Compatible results for the case of normal samples have been obtained by 
Williams [3], using another method. 
From the above results, the value of 


a 3(n® + 2n — 12)(n — 1)(n + 1) 
7 (n= 2)(n + 3)\(n+5) ’ 
is seen to approach normality as the sample size is increased. 

Inasmuch as the distribution of C for normal samples is limited in both direc- 
tions and is symmetrical, it is apparent that the Pearson Type II distribution 
may be considered representative. Fitting this curve to the moments given 
above, the equation of the frequency distribution is given by 


"ai 
y=w(1-S), 


_ (nt — n® — 13n? + 37n — 60) 
2(n? — 13n -+ 24) } 
* (n® + 2n — 12)(n — 2) 
(n? — lan +24)” 
r(2m + 2) 
a-2?™+1(T'(m + 1)} z 








m 


pe = 


The values of 62 for the distribution, for various values of n, are as follows: 


Sample size, n Bo 
5 2.300 
10 2.570 
15 2.684 
20 2.750 
25 2.793 


50 2.833 






Due to the effect of even moments higher than the fourth, the approximation 
afforded by the Type II curve is not reliable for samples containing less than 
about eight observations. As the sample size decreases below this limit, the 
extremes of the C distribution’ deviate increasingly from the extremes (a) 
of the fitted curve: with such a platykurtic distribution, therefore, the effect 
upon the lower significance levels vitiates the approximation. 

Although either 8, or the theoretical limits of the distribution of C could 
have been employed as a parameter of the fitted curve, it was considered ex- 
pedient to use the former. In any case, of course, the advantage to be gained 
would be in connection only with samples containing few observations (less 
than eight). The evidence afforded by empirical sampling indicates that use 
of the limits as a parameter might render the approximation less valid. 
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In order to facilitate use of the approximate distribution for samples of eight 
or more observations, the values of C associated with two probability levels are 
tabulated below in Table I. The ratio of each value of C to its standard error 
is also shown, to demonstrate the approach to normality. The significance 
levels recorded exclude 10% and 2% of the area under the curve, respectively. 
In most practical applications, these will be the 5% and 1% levels, respectively, 
since only positive values of C exceeding the tabulated value will ordinarily be 
considered significant. The tabulations were prepared from tables of the 
function I.(p, q) [5], where g = .5 and p = m + 1, with the transformation 

C 


zgz=l--— 


a?” 
TABLE I 
Significance levels of the absolute value of C 
Sample size,n P = .10 C.10/e¢ P = 02 C.02/e 

8 .0088 1.6486 .6686 2.1664 

9 .4878 1.6492 .6456 2.1826 

10 4689 1.6494 .6242 2.1958 
11 4517 1.6495 .6044 2.2068 
12 4362 1.6495 .5860 2.2161 
13 4221 1.6495 .5691 2.2241 
14 4092 1.6494 .5534 2.2310 
15 3973 1.6493 5389 2.2369 
16 .3864 1.6492 5254 2.2423 
17 3764 1.6492 .5128 2.2470 
18 .3670 1.6491 .5011 2.2513 
19 3583 1.6489 .4900 2.2550 
20 3902 1.6488 4797 2.2585 
21 3426 1.6488 .4700 2.2616 
22 3355 1.6486 4609 2.2647 
23 3288 1.6485 4521 2.2676 
24 3224 1.6484 .4440 2.2700 
25 3165 1.6484 4361 2.2717 
Normal (n = «) 1.6447 2.3262 


The distribution of C for normal] samples containing 20 or more observations 
is sufficiently normal, for most practical cases and for the more common signifi- 
cance levels, to permit use of a table of areas under the normal curve, in conjunc- 


tion with the standard error o, = \ / woah . Thed% significance levels 


shown in Table I result, at worst, in a one per cent error of probability estimate, 
if the normal approximation is used in their place: that is, if 1.6447 times the 
standard error is used instead of the tabulated significance level, the probability 
will be .0505 at most, for the values of n which are tabulated. 
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3. General discussion on the application of C. It may be wondered 
why the statistic C has been used, ‘rather than the more easily computed statistic 


n—l 
» (X; — Xess) 


Eat 
does not matter which is used, since C and C’ are linearly related. However, C 
may be regarded as symmetrically distributed about 0 in samples from a normal 
population to within at least four moments. Excessive departure of C from 0 
may be taken as indicative of the presence of non-randomness in the series, the 
actual significance test being based, of course, on the probability of obtaining a 
departure larger than a given observed one, under the assumption of a random 
series. Positive values of C, in general, correspond to positive correlation while 
negative values correspond to negative correlation between successive obser- 
vations. 

There are various ways of detecting non-randomness in a series of observations, 
such as regression methods, analysis of variance, etc. The use of regression 
methods implies that we must know in general the type of regression function 
to be tried. C is a very flexible statistic, on the other hand, for testing the null 
hypothesis that a series is random, no matter what the alternative hypothesis is. 
A thorough study of C as a statistic for testing the hypothesis of randomness in 
an ordered series should include a study of the power function of C for hypotheses 
specifying various types of non-randomness. However, we shall simply appeal 
to intuition in proposing the statistic C, and forego power function considerations 
in this note. In practice, the advantage of using C increases with the length of a 
series: lack of randomness in a single sequence of ten or less observations may 
ordinarily be detected by regression methods, in fitting a low order polynomial. 
In a longer sequence of measurements, on the other hand, the presence of com- 
plicated regression or of periodicity is often sufficiently obscured by variation 
to elude detection by any other than a flexible method. 

The statistic could be used to advantage in the field of applied statistics, in 
the investigation not only of variate series but of attribute series as well. For 
the latter purpose, an effort to tabulate the relationship between the level of 
significance and the percentage of either attribute would facilitate statistical 
investigation of random arrangement. A direct application could thus be made 
to binomially distributed attributes by a scalar assignment (0, 1) to the dichot- 
omy, followed by a procedure similar to that presented above. Similarly, the 
randomness of vectorial observations could be examined from the viewpoint of 
arrangement. The common method of treating such problems,—the ‘‘random 
walk method,””—has occasionally been found inadequate in dealing with specific 
forms of non-random order; this is especially true when the allocable cause of 
variation has a multi-directional effect. 

Needless to say, each of the fields of application considered so briefly above 
would require development before a routine, efficient method of investigating 
ordered arrangement could be established. Although probability level tables 


ce 


As far as a significance test is concerned, it clearly 
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have been provided in this paper for C as applied to normal samples, it is quite 
evident that tables for samples from other parent distributions would be needed 
for some of the applications mentioned above. 


4. An illustration of the use of C. Although one example has already 
been presented elsewhere [4] in which the distribution developed in Section 2 
has been employed, a typical application of the statistic to an example in the 
field of quality control will be given here in order to illustrate the mechanics of 
solution. The data presented in Table I] represent the percentages of defective 
product turned out daily, over a period of twenty-four days, by a single workman. 
The total output each day closely approximates five hundred parts: this fact is 
brought out to explain the calculation of x’ for the observed series of percentages, 
—it has no bearing upon the use of C. 


TABLE II 
Percentage of product rejected 

%, X X? 
7.4 54.76 
8.8 77.44 
11.4 129.96 
10.3 106.09 
11.9 141.61 
12.2 148.84 
10.0 100.00 
8.4 70.56 
9.4 88.36 
10.9 118.81 
9.9 98.01 
11.8 139.24 
10.0 100.00 
8.9 79.21 
9.7 94.09 
9.3 86.49 
12.0 144.00 
12.3 151.29 
10.3 106.09 
8.6 73.96 
10.4 108.16 
11.1 123.21 
9.4 88.38 
8.2 67.24 


Day 
1 
2 
3 
4 
5 
6 
7 
8 
9 


Totals J 2495.82 
nX° 2452.28 
te = 43.54 
C = .3636 (significant) x’ = 21.518 (23 degrees of freedom) (not significant). 
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The value of C derived from the data lies between the two significance levels 
tabulated in Table I; there is reason to believe that the data are ordered, or non- 
random. Computation of x*, however, has been carried out with the hypothesis 
that all product was made under thé same conditions (i.e. with a percentage 
defective equal to 10.108%, the mean of the group). The value so obtained is 
associated with a probability of about P = .50: the hypothesis is not disproved 
by this test. In short, the variability of the twenty-four observations could be 
considered random if it were not for the order of their arrangement. 
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ON CERTAIN LIKELIHOOD-RATIO TESTS ASSOCIATED WITH THE 
EXPONENTIAL DISTRIBUTION 


By Epwarp PAULSON 
Washington, D.C. 


Various likelihood-ratio tests and their distributions in samples from a popula- 
, “ Sadia 
tion having the elementary probability law —e Plt B< a < o, have been 
o 


studied by Neyman and Pearson [1] and Sukhatme [2]. In this note the power 
functions and the question of bias of several likelihood-ratio tests will be in- 
vestigated. The exponential distribution appears to be appropriate for dealing 
with problems involving the intervals of time between events which tend to be 
random, as for example the interval between consecutive telephone calls, or 
the interval between consecutive accidents to the same worker. 

To test the hypothesis H’ that the location parameter B is equal to some 
fixed value, it being assumed that the scale parameter o is known, we can for 
simplicity take the set 2 of admissible populations from which the sample might 
have been drawn to be {—« < B < + ©, o = 1}, while the subset w from 
which the sample must come when the hypothesis is true is {B = 0, o = 1}. 
Then the likelihood-ratio A, for testing this hypothesis is 


n 
- Zz % 
ee P(w max.) @ i=l on 
oe ren ae 
P(Q max.) eee 
C tml 


where x; is the smallest observation in a random sample of n. The region of 
acceptance of this hypothesis consists of all points in sample space for which 


Me < Ai < 1, 


1 
where X;, is chosen so that [ gi(A1) dy = 1 — a@, a being the level of significance 
Ate 
used and g (Ai) d\; being the distribution of \; when B is really equal to zero. 
The region \;. < Ai < 1 is equivalent to the region in the sample space for which 


0<n<hjk = ~ Me. 


For any value of B the distribution of x; is known [8] to be 


$1 (21) dx, = ne n(z,—B) dz, : 
301 
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poo ares 
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Setting B = 0, the relationship between k, and a is 
ky 
| ne“ dzj=1l—a, so e™“'=a@, 
0 
When B < 0, the power function P(B), for this test is 


Ky 
P(B) =1—- | ne" dz, = 1 — e™[1 — al. 
0 


ky 
When 0 < B<k,, P(B) =1- i ne") dz, = ae”. When B > hk, 


B 
P(B) = 1. 

Since e’” > 1 if B > 0 and also e™ < 1 if B < 0, P(B) is obviously > a if 
B #0. This test is therefore completely unbiased in the sense of Daly [4]. 
In addition, it is not difficult to prove that this test has the unusual property 
of being a uniformly most powerful test with respect to all alternatives. 

To test the hypothesis H” that the location parameter is equal to some fixed 
value, say B = 0, when the scale parameter o is unknown, the likelihood-ratio 
is easily seen to be 

Zz (a; — 21) 


Ne = | —— ta cea 
> 2; a 
t==1 
a (x; — %) 


The region of acceptance consists of all points in the sample space for which 


n 











1 
Nee < Ae < 1 where [ go(A2) dXg = 1 — a. This is equivalent to the region 
Age 


(1) 0< oo Ske; k, = (n — 1) (1 — dae") <a") 
Ze (x:—21) Age 


The relation between kz and a is easily found from the distribution of t when 
B = 0, which is known to be [3] 


at 
. 
[r+ | 


ke k —(n-1) 
Therefore [ g(t) dt = 1 — a, so [1 + z : | =a. 
» - 


It is somewhat easier to find the power function of this test by considering the 
region of acceptance as made up of points in the zx, , s plane for which 


. (x5 — 21) 


tml 


n—l 


go(t) dt = 








kes 


osa35— where s = ; 


which is identical with the region in (1). 





LIKELIHOOD RATIO TESTS 


The joint distribution of z; and s is [3] 


¥ilti, s) dx, ds = 3(x1) dxy-d,(s) ds, 


where 


_ ‘aid 
$3(21) da, = — e "1" da, 
a 


ie n—1 
(" *) gs” 2 ele ds 
ga(s) ds = = 


(n — 2)! 
When B < 0, the power function P(B) of this test is 


o koa/n 
P(B) =1- [ de I pind) dae i - "A - a 


When B > 0, the power function is 


P(B) =1- | - as [ “_— 
(2) Bn/ko B 


on aenP'* 4 1 as Meee | = ae" n a 1,24" = +m), 


where I[p; z] = 


which is the form in which the Incomplete Gamma Function has been tabulated 
(5). 

Since o must be positive, e"’’” < 1 if B < 0 and therefore P(B) > a in the 
interval —«» < B <0. Toshow that P(B) is > ain the interval0 < B < o, 
it is simpler to work with the expression for P(B) as a double integral in (2), 
than to differentiate the power function directly. Performing the integration 
with respect to 7, 


P(B)=1+) [ee " _ 1). g4(s) ds. 
nlke 
Differentiating with respect to B, 
P’(B) = [ . = eel (5) ds. 
nike 


The integral expression for P’(B) is obviously positive. Therefore since for 
B > 0 the derivative is always positive the function must be monotomically 
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increasing in this interval (0 < B < +),so P(B)is > awhenB>O. There 
fore this test is also completely unbiased. 

We now consider the hypothesis H’” that two samples are drawn from ex- 
ponential distributions with the same location parameter, assuming it is known 
the samples must have come from two exponential distributions with the same 


. Eten 
scale parameter. Given a sample of m; values of x drawn from — e 77”! dz 
o 


and another independent sample of nz values of y drawn from : e ¥-80/* dy the 


hypothesis we wish to test is that B, = B,. Let x be the smallest of the n, 
values of x and y; be the smallest of the nz values of y, let L be the smallest of 
the nm + nm, = N values of both xz and y. Then the likelihood ratio for this 
hypothesis is 
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where z= nly — 24), if y> 2% 


= ny (2, = Yi), if T1 > Y1, 
n2 


and u= » (x; — x1) + > (yi = y1). 


t=l1 


The region of acceptance, As. < A3 < 1, is equivalent to the region 0 < Z < Kgu, 
where K; is again a function of a, the level of significance, the exact relation 
being 
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It is known [8] that u is independent of Z, and that its distribution is 
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The distribution of z is somewhat complicated; but it can be derived by observ- 
ing that the probability that z lies in any infinitesimal interval z: + } dz is 
the sum of the probabilities that ne(y: — 21) and m(z, — y;) lie in that interval 
and by then using standard methods for finding the distribution of the difference 
of two variates. For the case G = B, — B, > 0, the distribution f(z) of z is 
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For the case G < 0, the distribution of z can be derived from (3) by interchang- 
ing m and nz and putting —G in place of G. 

The power function of this test can now be derived. For the case G > 0, 
the power function P(G) is 


2 noG noGlk3 noG 
P(G) =1 -{ I du | fule)os(u) dz — I du [ fulz)oo(u) de 
(4) < kgu 
+ [ ange dg fe) ostu) ic} 


Upon integrating out and simplifying, the power function becomes 


—n,Gie ‘ 
re =o) + rv 2 2 
kso 


Nm + Ne 


ae) e l mG(1 + | 
+0(™ ay 1-—I|N —— 


Ne pon Ne y | =~ % G(n, — mks) 


NN + Ne Nn. — mks kso 


The power function when G < 0 is easily derived from that for G > 0 by every- 
where interchanging n; and ne and substituting —G for G. 

To show that P(G) > a when G = 0, it is only necessary to show that the 
derivative P’(G) of the power function is always positive when G > 0, and al- 
ways negative when G < 0. It is again considerably simpler to use the expres- 
sion for P(G) as a double integral. For the case G > 0, integrating with respect 
to z in (4), 











PG) = 1 — —2— fa — oir] 
Mm + Ne 


noGlk3 me 4" ' \ ‘ 

nyz a n 

+ [ aime) Oo etude 
1 2 


(me + me") | sie 
— —<-—" o 
2 3 


where [f(z)]? = f(b) — f(a). Upon differentiating and simplifying, 
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Both integrals are easily seen to always be positive, so P’(G) is positive when 
G>0O. Inthe same manner it can be shown that P’(G) is negative when G < 0. 
Therefore this test is also completely unbiased. 
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The question of investigating the bias of the likelihood-ratio tests for (a) 
testing the hypothesis that ¢ = oo when B is known and (b) testing the hy- 
pothesis that ¢ = oo , nothing being known about the value of B, are practically 
identical with the analogous problems for a normal distribution. The results 


are also the same, for the \ test for (a) is completely unbiased, while that for 
(b) is biased. 
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1. Introduction. The number of mathematically significant figures in the 
solution of simultaneous linear equations has received attention from a number 
of writers [1-6]. It is an important subject, not only in least squares and 
correlations, but in many other problems of science where simultaneous equa- 
tions arise: it may not be amiss, therefore, to examine it from a fresh start, 
particularly since (as will be shown) some of the rules that have been published 
on it fail in certain frequently occurring circumstances. 


2. Definitions. Before proceeding into the subject it will be necessary to dis- 
tinguish between the computer’s terms “significant figures’ and “determinate 
significant figures.”” The former are the figures that compose a number, without 
the consecutive ciphers that precede or follow them, merely to locate the decimal 
point. ‘‘Determinate significant figures,” on the other hand, are figures that 
are justifiable on computational grounds. From the computer’s point of view, 
the number of significant figures remains independent of what is statistically 
significant. To avoid confusion in what follows, the term “significant figures” 
will be used in the computer’s sense, and the adjective ‘‘determinate”’ will be 
supplied where mathematical determinacy is implied. 

To avoid prolixity the term ‘observational error’’ will include any uncertainty 
arising either from errors in the observations or from the statistical nature of 
the problem (e.g. sampling errors, grouping errors, etc.). The observational error 
of the result is independent of the particular sequence of computation followed and 
the accuracy with which it is carried out. 

The term “computational error” will include all the additional uncertainties 
arising from the approximations occurring in the particular sequence of computa- 
tion used, including the ‘‘rounding off” of the final result. The computational 
errors, unlike the observational errors, depend in general upon the sequence of the 
intermediate steps used in the computation as well as on the number of significant 
figures to which they are carried. 








3. Criterion of an adequate computation. If the number written down at the 
end of a computation is to serve its purpose the maximum possible computational 
error must be suitably limited. 


A decimal representation of a number containing f significant figures is subject 
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to an uncertainty (upper limit of absolute error) of 5 in the (f + 1)th place. 
It has, therefore, a possible relative (not absolute) error of representation some- 
where between 5 X 10°”*” and 5 X 10°’, in magnitude. This relative compu- 
tational error sets the limit to any valid final rounding off. Regardless of the 
accuracy to which the intermediate steps of the computation have been carried, 
this relative computational error introduced by the final rounding off alone 
must be suitably limited. 

In case all of the accuracy obtainable from the data is not needed in the result, 
the sum of the maximum possible computational error (including the error of 
the final rounding off) and the maximum possible observational error must be 
kept below the error which can be tolerated in the result. 

In case all of the accuracy obtainable from the data is needed in the result, 
the maximum possible computational error in the result (including the error of 
the final rounding off) must be negligible in comparison with the uncertainty 
(observational error) in the result arising from uncertainty in the data. Just 
how small a fraction of the observational error 1s ‘‘negligible”’ 1s necessarily a matter 
of judgment, and will depend upon the nature of the problem. A computational 
error that would be wholly negligible in some ordinary computations might be 
intolerably large in the adjustment of an accurate geodetic survey. In any case 
the only basis for a valid judgment of the adequacy of the computation lies in a 
comparison of (i) the maximum possible computational error that can arise in 
the sequence of computations including the final ‘‘rounding off,”’ with (ii) the 
observational error of the result arising from the observational errors inherent 
in the data. 
















4. Propagation of error in a system of linear equations. Assume that 


(1) z Att: = bs, s=1,2,-+-,n, 
t 
















is a set of simultaneous linear equations derived in some way from observations 
and in which the coefficients a,, and the absolute terms b, may all be subject to 
observational error. If the relative (not absolute) observational error of a 
quantity q be represented by 6, it may readily be seen that 


6x; = he 2X (x4. /2;)AnjOnx 50n% + ke (b,/2;)Asj5b, 


(2) 
X 2X A nk One 5Onk 


6A 


1 
where A is the determinant of the coefficients a,,, and Ay, is the term corre- 
sponding to a, in the reciprocal (not the adjoint) determinant. 


5. Upper limits to observational errors. The sign and magnitude of the 
relative errors 6a,, and 5b, are unknown, but we shall assume that it is possible 




















SIGNIFICANT FIGURES 


in any problem to assign to them upper limits 


| dan. | and | 6b, | 


which in magnitude they cannot exceed. If the problem is such that the values 
of each of the da, and the 6b, are wholly independent of each other, it is then 
possible that their magnitudes may all reach their upper limits | 5a), | and | 5b, | 
simultaneously, in which case upper bounds of 6x; and 6A may be placed at 


| 6x;| = a » | (wi./2j)Anjane | | Sane| + Do | (b./23)Aaj| | 6. | 


3 
” }5A| = 22 De | Anwane| | das | 


6. Indefiniteness of the problem in the general case. Tie values of the dan 
and 6b, may not be independent of each other, in which circumstance knowledge 
of the law of their dependence would make it possible to assign upper limits to 
the magnitudes of 6z; and 6A. These upper limits can not be larger than the 
upper bounds shown in equation (3), and in special cases they will be much 
smaller. Since the dependence of 6a,, and 6b, may in general have any form 
whatever, cases can and will occur in which the upper limits of the relative 
errors of 6x; and 6A may have any ratio whatever. 


7. Case of independent errors. Any general discussion of the errors that can 
occur in x; and A must be based either on some special assumption or on the 
limiting assumption that the errors are independent. It is this latter assump- 
tion that underlies the usual discussion, and will be the basis of what follows. 
Equation (3) gives the upper limit to the 6x; and 6A under these assumptions. 


8. The ratios of | 5z;| and | 5A| are still indefinite in spite of the assumption 
of independent errors in the coefficients. However, equation (3) does not deter- 
mine any definite ratio or inequality between the upper bounds | 62; | and | 6A |. 
The nature of the observations may be such that some of the errors in the an 
and b, are very small and some relatively large. Not infrequently it is safe to 
assume that some of them are free from appreciable error and to ascribe all the 
error of the x; to the error in one or two of the a,x, orb,. If any statement of a 
definite relationship, either as an equality or an inequality between | 6A | and 
the | 6x; | is valid for all possible sets of linear equations, it must at least hold 
in the special case in which the errors of all the b, and the errors of all except one 
of the ay, are negligible. 

If such a statement of a definite general relationship between these upper 
limits of errors can be made, it must be possible to write down an equation or an 
inequality between any one of the expressions | Ay. | and some or all of the 
corresponding expressions | (x;/2;)An;|,j = 1, 2, ---,, that will remain true 
no matter what be the values of the a,, and the b, in the original set of simulta- 
neous equations. It is obvious that the ratio of | An| and | (x:/x;)An;|, 
(j * k), depends upon the values of the aa; , and sets of equations can be found 
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to give any assigned value to that ratio. It is therefore impossible to state any 
rule that will restrict the ratio of the relative error of A and the relative error 
of any one of the z;, valid for all possible sets of linear equations. 


9. Definite statement about the sum of the relative errors in the unknowns. 
However, in the summation A | dx; | there occurs the term corresponding to 


2 
j = k, for which | (x,/x;)An; | = | Anz |, so that under the assumption that the 
ay, and b, are independent sources of error, we may write the inequality 


(4) dD | 6x;| £ |5A| 


i 


which states that the sum of the upper bounds to the relative errors of all the z; 
cannot be less than the upper bound to the relative error of the determinant A, 
A corresponding statement can easily be proved for the standard deviations. 

A limiting case can be constructed in which the inequality (4) reduces to 


(5) 2, | 6x;| = | 84 


and in which all of the | 6x; | are equal. For this case, 
(6) |54 | = n| 62; | for all values of j. 


If n < 10 it is obvious that there will be at least one more determinate signifi- 
cant figure in each of the z; than in the determinant A of the coefficients. 

It is frequently assumed that the number of determinate significant figures in 
the solution for any unknown cannot exceed the number of determinate signifi- 
cant figures in the determinant A of the coefficients. We see now that this state- 
ment can not be generally valid, even under the assumption that the a, and b, 
are independent sources of error. As a matter of fact, it is necessary in some 
cases to compute some or even all of the unknowns to more significant figures 
than are determinate in the determinant A of the coefficients, if one would retain 
in the result all the accuracy that is obtainable from the data. 

Cases in which the relative observational error of every one of the unknowns 
is less than the relative error of the determinant A probably occur rarely in 
practice; in fact the only ones that I have seen are those that I constructed 
purposely to show that such a thing is possible. However, cases in which the 
relative errors of one or several but not all of the unknowns are much smaller 
than the relative error of the determinant A, occur fairly frequently. 


10. Remarks on the case of “near indeterminacy.”. The major interest in 
curve fitting centers around the condition of ‘near indeterminacy,” i.e., of a 
small or near vanishing determinant A. Even in the circumstance where the 
relative error of the determinant is much greater than the relative error of some 
or all of the coefficients and absolute terms, the relative error of one or more of 
the unknowns may be much smaller than the relative error of the determinant, 
as may be seen from what follows. 
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In accurate experimentation the endeavor is, wherever possible, to arrange the 
experiment so that the quantity sought comes directly from the measurement as 
represented by an equation such as 


(7) x=D. 


However, so ideal an experimental arrangement is rarely if ever possible, and it 
is a common experience to find that the measurements are represented by an 
equation such as 


(8) rct+qtrze+sut+.--- =p, 


where gy, 7z, su, etc., are small corrections that must somehow be evaluated. 
For simplicity, the discussion will be confined to the almost trivial case 


(9) x + qy = Pp. 


Not infrequently the only way the correction can be evaluated is to rearrange the 
conditions of the experiment so that another equation is obtained in the form 


(10) zr+qy=p’. 


Sometimes the nature of the experiment is such that it is not possible to change 
the coefficient of y by more than a small amount, under which conditions 


(11) q =q(1 +8), 
and 
(12) . p’ = p(l + a), 


where 6 and @ are small in comparison with 1. The solution of equations (9) 
and (10) now gives 


py’ — pg 
13 : = p(1 — : 
(18) v7 p(1 — «/B) 


The quantity g’ — gq seen in the denominator of this equation is the determinant 
A of the coefficients, and by equation (11) its value is Bg. Since 8q is assumed 
to be small here, the solution for x encounters a near vanishing denominator. 
It would, however, be wrong to assume that the number of determinate signifi- 
cant figures in x that can be obtained by solving the equations is necessarily 
limited to the number of determinate significant figures in the denominator A. 

If the experimenter has been fortunate in finding suitable experimental condi- 
tions, the denominator A = §gq, although small in comparison with either q’ or q, 
will still not cause difficulty. It will be observed that the coefficients of q’ and q 
in the denominator are equal (both being unity). Now if the coefficients p 
and p’ in the numerator are nearly enough equal, so that g’ and g occur in both 
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numerator and denominator so nearly proportionally that the uncertainties in 
q and q’ produce nearly compensating errors in both numerator and denominator, 
then zx will be given to more determinate significant figures than are found in 
the denominator A. It can then be said that the experiment is successful in 
evaluating the correction term gy in equation (9). 

On the other hand, in less fortunate circumstances, to the exasperation of the 
experimenter, the denominator A = q’ — g = gq is not only small, but p’ and p, 
although still nearly equal, differ enough so that the errors in q’ and q are not 
compensated by the nearly equal coefficients in the numerator. The experiment 
will then fail to improve the approximation p for z by failing to evaluate the 
small correction gy in equation (9). This would be an inherent defect in the 
experiment and could not be removed by any manner of computation. 

The same conclusion would of course be drawn from the coefficient of p (viz., 
1 — a/®) at the extreme right of equation (13). It is not the size of 6 that 
alone determines the number of determinate significant figures in z, it is rather 
the ratio between a and 8. In the fortunate experimental circumstances de- 
scribed above, the near equality of p’ and p offsets the near equality of q’ and q¢ 
by reducing the term a/8 to a value small compared with unity; the term a/8, 
being small, acts to reduce the effect of the uncertainties in q and q’ (i.e., in q 
and £) in the evaluation of x. On the other hand, in less fortunate circum- 
stances, the correction term a/8 can not now shield x from the uncertainties in ¢ 
and q’ since the relative difference a between p and p’ is not small enough to 
reduce a/8 to innocuity. 


11. Numerical illustration of compensating errors. As a “horrible example” 
especially constructed to emphasize the theoretical possibilities, take the fol- 
lowing special case— 


1000.10000z + 10.00000y = 1010.10000 
1000.00000z + 10.00000y = 1010.00000 


(14) 


wherein it is assumed that the coefficients and the absolute terms (assumed to 
be derived from the observational data) are all correct to the fifth decimal place 
as given, and no closer estimate of their errors is possible. So far as known, the 
upper limit to the absolute observational error of each is then the same, i.e. 
5 X 10°, but the coefficients of x (ay, and a»), and the absolute terms (b; and be), 
all have nine determinate significant figures, while the coefficients of y (an 
and dz), have only seven. Thus, 


| 6a | > 5 X 10°, | dan, | > 5 X 10°, |db1| 5 X 10°”, 
| dbo| > 5 X 10°, 


| day | > 5 XK 10”, | daze | > 5 X 10°”, 
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and x = 1, y = 1, A = 1, whereupon a substitution of values from (15) into (3) 
gives the inequalities 


(16) |ér|>3 xX 10%, |dy| +3 X 10°, |54| > 1.01 x 10”. 


So far as known, the determinant A may thus be in error by as much as 1 per 
cent, and y by as much as 3 per cent, vet x is known closer than 1/30th per cent. 
Here the value of the unknown 2 cannot be adequately represented by less than 
four significant figures, and might even require five, in spite of the fact that 
neither A nor y requires more than three significant figures to represent all that 
is certainly known about them. 

The reason for this disparity in relative errors can be more easily seen by 
substituting numerical values for all the coefficients in the expression for x 
except dic and dy. The possible relative errors of aj and de are, as noted 
above, about 100 times as great as the possible relative errors of ay, da, bi, 
and b. , and are the controlling errors in A. In the solution 







x = 1010.10000a2 — 1010.00000a12 
~~ 1000.10000a22 — 1000.00000a12’ 


however, both a2 and de occur in both numerator and denominator, and more- 
over the coefficient of each in the numerator is nearly equal to its coefficient in 
the denominator, so that a change in either aj. or a2. changes both numerator 
and denominator nearly proportionally, with the result that their ratio z is 
known much more accurately than either the numerator or the denominator A. 

This kind of compensation of errors in a computation is not confined to the 
solution of simultaneous equations (and it is not an infrequent occurrence in 
other computations). This is one of the many reasons why it is impossible to 
give general rules for the retention of significant figures that will be valid for 
all types of computations. 


12. Geometrical analogy. Moulton [4] illustrated his reasoning by the fol- 
lowing geometrical analogy. The solution of three linear equations is equivalent 
to finding the point of intersection of three planes. When the determinant of 
the coefficients is small in comparison with the coefficients themselves, these 
planes are either nearly parallel, or the line of intersection of any two of them 
is nearly parallel to the third. In these cases small uncertainties in the location 
of any one of the planes correspond to large uncertainties in the position of their 
point of intersection. 

In the first circumstance the planes might all be nearly parallel to one of the 
three coordinate planes, with the result that large uncertainty would afflict the 
value of the determinant and two of the unknowns, the third being much more 
accurately determined. 

In the second circumstance, the line of intersection of two of the planes might 
be nearly parallel to one of the coordinate axes. When that happens, large un- 
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certainty will afflict the value of the determinant, but only one of the unknowns, 
the other two being much more accurately determined. 

This geometrical analogy can be extended to cover simultaneous equations 
with any number of unknowns. Near-vanishing of the determinant A of the 
coefficients necessarily implies relatively large uncertainties in the determinant 
and also in at least one of the unknowns, but not necessarily in all of them. 
These are, of course, very special cases, but, as noted above, they are of frequent 
occurrence in actual problems. 


13. Evaluation of computational error. The relative computational error in 
x; must be kept within certain definite limits which depend upon the particular 
problem to be solved (section 3). To do this it is necessary to be able to calcu- 
late an upper bound to the relative computational error inherent in any particular 
sequence of computations. 

In many computations it is easy to write down a simple formula that will set 
an upper bound to the relative computational error involved in that particular 
sequence. This formula contains numbers f; , fo , fs , etc., each representing the 
number of significant figures accurately computed at some particular step. 
Once a simple formula for relative computational error is written down, it is 
easy to choose values of fi, fo, fs, etc. that will give an upper bound to the 
relative computational error not larger than the permissible limit of maximum 
possible computational error outlined in section 3. This method of determining 
an upper bound of the relative computational error should be used whenever such 
a simpl2 formula can be found. For example, to compute x from equation (13) 
we may use the following sequence: 7; = q' — q, re = 11/q = 8,73 = p’ — py s = 
13/p = a, 1s = 14/72 = a/B, 76 = 1 — r5 = 1 — @/B, 77 = pre = (1 — a@/B) = 2. 
xz may then be written as a function of these partial results, viz.: 

(18) z= 1, = pre = p(l — 7s) = p(l — 14/72) = p(l — 13/pre) 
= p(l — mg/n). 
Applying first order error theory we find 
_a/6_ | 


le(x)| < —_ 


(19) i— a/b | {| e(ry) | + | e(re) | + | e(rs) | + | e(rra) | + | (rs) |} 


+ | (76) | + | €(r2)| 


where e¢(r;) represents the relative error in r; arising from the computation by 
which r; was determined from the precéding partial results, 71, 72, --- ,7i-1, 
and ¢(z) is the total relative computational error in x when so computed. It 
is easy to keep e(x) within any desired limits by suitably limiting each error term 
of (19). Since a computation accurate to f significant figures involves a relative 
computational error not greater than 5 X 10°’, any desired limits can then be 
set to each error term of (19) by a proper choice of the number of significant 
figures that should be carried in that step. 
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Unfortunately there seem to be no reasonably simple formulae for determining 
upper bounds of the relative computational errors that arise in the solution of 
simultaneous linear equations in more than two variables. This does not ab- 
solve the computer from the necessity of ensuring that his computational errors 
are suitably limited. 

The method I have found most economical is to carry the solution of simulta- 


neous linear equations to the capacity of the machine, and as each partial result 
r; is obtained, write it as 


ri(1 + €i), 


where 7; is the value actually found and ¢; is a positive number representing the 
accumulation of uncertainty introduced by all preceding steps in the computa- 
tion. At the end of the computation each of the unknowns is found. in the form 


(20) x(1 + €)), 


where x; represents the value found and e; is the upper bound of the relative 
computational error in z;. 

A comparison of e; with the upper bound of the observational error | 6z; | of 
equation (3) will then indicate whether the computation is adequate. If the 
comparison shows that the computation was inadequate, it will show in which 
steps the number of significant figures f; was too small, and by how much. 
The computer can recompute, carrying these steps to the requisite number of 
figures with the assurance that his recomputation will then be adequate. The 
comparison will further indicate in which steps if any the number of significant 
figures f; was larger than necessary. 

When a computer has thus set suitable upper bounds to the relative computa- 
tional error in the solution of a set of linear equations, he is in a position to plan 
solutions of future similar sets so as to perform his computations more eco- 
nomically and yet safely. This is especially true when the solution of simulta- 
neous linear equations arises week after week in routine testing. 


14. Conclusions. Summary rules have been published, purporting to be safe 
guides to computers in avoiding needless work, and ensuring that the computa- 
tions are carried to a sufficient degree of accuracy. Many of them are useful 
guides for certain types of computation and for limited ranges of the numerical 
values entering into the computation, but none of those that I have seen can be 
used generally. The only safe rule, where the matter is of importance, is to 
calculate the maximum possible computational error that can enter in the par- 
ticular sequence of computation followed, and make sure that it is kept within 
the necessary limits. 

It is sometimes necessary to carry the intermediate steps of a computation to 
many significant figures beyond the significant figures given in the data, or kept 
in the result. The relative error of one of the unknowns may be very much 
smaller than the relative errors of the data from which it is computed, while the 
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relative error of another of the unknowns may be larger. The methods of 
ensuring that the computations are adequate are outlined in section 13. 

For the best sequence to follow in the elimination of the unknowns, I shall 
pass along a suggestion of Dr. W. Edwards Deming which he gave in one of our 
discussions of this subject. I venture to pass it along, because it has worked in 
every special case that I have constructed in an attempt to prove that it does 
not hold generally. If ever the suggestion fails, the computer may change the 
sequence; but in any case he is obliged, as stated above, to calculate the maximum ° 
possible computational error that can enter into his calculations. Dr. Deming’s 
suggestion is this: ““To evaluate some but not all of the unknowns to the highest 
possible computational accuracy, retaining as few significant figures as possible 
in the intermediate steps, solve the equations by successive elimination, elimi- 
nating first and evaluating last the unknowns of greatest inherent relative 
accuracy.” 


15. Summary. Expressions are given for the maximum observational error 
in the unknowns of a system of simultaneous linear equations, in terms of the 
relative errors of the coefficients and absolute terms therein. In order to extract 
all the information possible from a system of linear equations representing ob- 
servational results, it is not sufficient in general to assume that the relative errors 
in the unknowns are as large as the relative error in the determinant of the 
system. In many problems the computation of some of the unknowns must 
therefore be carried to more significant figures than are determinate in the 
determinant of the system. Methods are outlined for evaluating computational 
error in the solution of linear equations to ensure that the computations are 
adequate. 

In conclusion I wish to express my thanks to Dr. W. Edwards Deming who 
has given much of his time to assist me in the preparation of this paper. He has 
made valuable suggestions on the material to be included and the general manner 
of presentation. In addition he has criticized the manuscript in detail and 
assisted in the final revision. 
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ON MECHANICAL TABULATION OF POLYNOMIALS 
By J. C. McPHERSON 


International Business Machines Corporation 


1. Introduction. The purpose of this paper is to show how automatic 
accounting machines, which have been used previously in evaluating such 
quantities as =x" and =x" "y, may be used in the preparation of mathematical 
tables of integral powers, of polynomials, and of functions which can be approxi- 
mated by polynomials. These tables may be prepared for any desired intervals 
of the argument such as 1, +5, rds, 4, 4, ete. 

The method is an adaptation of the general theory of ‘cumulative’ or “pro- 
gressive” totals which has proved useful in computing moments and product 
moments both with and without accounting machines. The reader unfamiliar, 
with the mathematical method and its machine applications might refer to such 
presentations as those of Hardy [1], Mendenhall and Warren [2, 3], Razram and 
Wagner [4], Brandt [5], and Dwyer [6, 7]. The main feature of the method is 
the computation of summed products or of summed powers by means of succes- 
sive cumulated additions. It is shown in this paper how it is possible to use 
this same process in constructing tables of powers and tables of polynomials. 


2. The Cumulative Formulas. If the numbers F; are defined and finite 
forz = 1, 2,3, --- , (a — 1), a, and if these values of F, are cumulated for z = 
a,x = a — 1, etc., then the value in the row headed by x = 1 can be written 
as 'T,. If these cumulations are cumulated successively with the superscript 
indicating the order of the cumulation and the subscript indicating the value of x 
which heads the row, then 


9 


“T; = 22F,, a = z a Fee a = Z ea 1) F., 


ip, = pet 2+ Ve 


3! Me 


and in general for 7 < j, 


ip — vit $9 —- G+ 019 
(1) Py = ay 


Formula (1) is basic to much of the previous work involving cumulative totals. 
Various authors have studied such important special cases as (A) where F; 
equals the frequency function f, , (B) where F, = xf, , and (C) where F, equals 
the sum of all the values of y having the same xz value. These special cases have 
been found very useful in computing moments and product moments. 
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The moments may be expressed in terms of the cumulations in a variety of 
ways. The diagonal formulas have the differences of zero as coefficients and are 
expressed in terms of ‘T;, “71, “T2, “T:, °Ts, etc. The columnar formulas, 
whose coefficients have been recently studied [6, 7], are expressed in terms of 
cumulations of the same order, ’7,, with j fixed. Razram and Wagner [4] 
have given formulas which utilize the entries of different rows and different 
columns but which demand fewer entries for the formulas. Razram and Wag- 
ner worked out the formulas through =z‘f, but the argument holds for Dz‘F,. 
For purposes of comparison the values of =z'F,, i = 0, 1, 2, 3, 4, as they appear 
in the diagonal, columnar, and Razram-Wagner systems are presented in 
Table I. 
TABLE I 

Values of 2x'F, fori = 0, 1, 2, 3, 4. 





Fz Diagonal 





Columnar Razram-Wagner 


=F, oY 17; iT’, 

22F, 27; 27; 27, 

r2*F, | 27, + BT 2 7, + 7's 37, ++ ATs = FT igs 
=2°F, 27, oa 6°T'» a 617’; 4T; “Bb 4“‘T, + ‘T’; 27; + 67’. 


Sr'F, | 271 + 14°72 + 3647's + 24574 | 8T) + 11572 + 1157's + 874 PT i420 + 1257243 


In developing the theory of the later sections of this paper I have developed 
further formulas of the type shown by Razram and Wagner since these formulas 
have fewer terms than do those of the other systems and the coefficients are 
factorable by (j — 1)!/2. These formulas for 22°F, , with s even, feature such terms 
as *T, + °T2 = *Tise, Tei, etc., so that there are two entries from the same 
column. For the purposes of this paper it is preferable to have a single entry 
from each column and this situation results from continued application of the 
formula 


(2) Tony = TT: + Tan = OT + 2 Tin. 


The formulas for 2z*F, with s S$ 12 are given. The alternative forms are given 
for the formulas involving even values of s. 


=F, =°T,, 2aF, = °T), Da *F, = *T, + Te = *Tigg = *T1 + 2°72, 
2x i = 7, + 6 'T,, 2x 7. = Vass + 12 "Fons 
= °T, + 2°T, + 12‘T. + 24 Ts, 


Sa °F, = °T, + 30 ‘T, + 120°73, 
rr °F, = *Ti2 + 60 “Tors + 360 "Ts44 
= °T, + 2°7. + 60 ‘T. + 120 °7s + 360 °T: +- 720 "7, 
ta 'F, = °T, + 126 ‘T, + 1680 °7; + 5040 *7,, 
Da *F, = "Tie + 252 *Tei3 + 5040 "Tas + 20160 °T 145 


= "7, + 2°7, + 252 ‘7, + 504 7; + 5040 °7; + 10080 "Ty 
+ 20160 *7, + 40320 °7;, 
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°T, + 510 ‘T, + 17640 °T; + 151200 °7, + 362880 “T;, 
= *T 142+ 1020 1243 + 52920 'T's44 + 604800 °T44s + 1814400 “Ts 
°T, + 2°T, + 1020 ‘7, + 2040 *7; + 52920 °T; + 105840 ‘7, 
+ 604800 *7', + 1209600 °7’, + 1814400 7, + 3628800 "7, 
°T, + 2046 ‘T, + 168960 °7; + 3160080 *7, + 19958400 “7; 
+ 39916800 "7, 
®T 142 + 4092 °T'43 + 506880 77314 + 12640320 Tass 
+ 99792000 “T's1s + 239500800 “7's,7 
°"T, + 2°T, + 4092 *T. + 8184 °7's + 506880 °7'; + 1013760 "7, 
+ 12640320 °T, + 25280640 °7, + 99792000 “7; 
+ 199584000 "7, + 239500800 “7, + 479001600 “7;. 
The derivation of these formulas is obtained with the use of (1), with the use of 
(4) T= Tin + T:, 
and with the use of formulas of lower order. For example we have from (1) 


y (x + 4)(x + 3)(x + 2)(x + I)a F. 


at 
120 or 


so that 
Sax °F, = 120°T; — 10 Sr ‘F, — 35 Dx °F, — 50 Dx ’F, — 24 UF, 


which after substitution of Zr‘F,, 22°F, , etc. and simplification results in the 
value ?7, + 30‘72 + 120°73. 


3. Tables of powers. If F. = 1 when x = a, but is zero otherwise 
then =z ‘F, is equal to a’. It follows that the value of a* can be obtained from 
the successive cumulations of this F, with the use of (3). For example in 
Table II 


TABLE II 
Cumulations of F; = 1, when z = 6, 
0, when x ¥ 6. 


ee 


ONMaA RON & 
cooooocor 
ONDER whe 
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6° 


°T, + 2°T, = 6 + 2(15) = 36, 
6° = °7, + 6*T. = 6 + 6(35) = 216, 
6 = 77, + 2°77. + 12°7, + 24°73 = 6 + 2(15) + 12(35) + 24(35) = 1296, 


The values of °7; , °7'2 , *7'2 and °7 for a = 6 are italicized in Table II. 

To get the values of 5’, 5°, 5*, etc. it would be necessary to start to cumulate 
from x = 5. Now since the values of ‘7; are unity, it follows that the values 
for a = 5 can be found by taking the entries above those for a = 6. Thus 
°T, = 5, *T2 = 10, ‘T, = 20, °T; = 15 with 5 = 5 + 2(10), 5° = 5 + 6(20), 
5* = 5 + 2(10) + 12(20) + 24(15). It is evident in general that the values for 
any a’, a’, a’ can be obtained by taking the row headed by a as the bottom row. 
Thus using a = 8, we have 8° = 8 + 2(28), 8° = 8 + 6(84), etc. It then appears 
that we may omit the z column of Table II and consider the cumulations to be 
ascending cumulations for a rather than descending cumulations for z. 

A more satisfactory course is to cumulate the coefficients so as to eliminate 
the multiplications. Thus the value of 6’7; could be obtained without multi- 
plication by cumulating 6, 0, 0, 0, 0 --- rather than 1, 0, 0, 0, --- . Several 
cumulations may be carried on at the same time so that the additions are not 
necessary and the tabulation results in a table of the desired powers. 

In preparation of a power table, the formulas (3) become a series of instruc- 


tions on the way in which we are to do the cumulating. For instance the 
fotmula: 


x’ = 5040°7; + 1680°7; + 126°T, + *7,, 
tetis us that to form a table of the seventh power we must cumulate’ the coeffi- 
cient 5040 eight times; add in the coefficient 1680 when there are six operations; 
the coefficient 126 when there are four; and the coefficient 1 when there are two 
remaining. A change in subscript tells us that the coefficient when first included 
forms a separate total ahead of the ones already partly figured. When the sub- 
script does not change, the coefficient is to be included in the first summary card 
total. The final cumulating operation prints the actual table. 

To prepare a power table by machine we secure a set of cards punched all alike 
with the numbers from 1 to 9 punched diagonally in successive columns across 
the card. The machine is wired to add the coefficient of the highest term by 
selecting the proper digits from the diagonals, cumulate after each card and sum- 
mary punch each total. This way of starting saves one cumulation. The 
summary cards are cumulated repeatedly in the same manner until the number 
of operations indicated by the highest term is completed. When the number of 
operations remaining equals j of another term ’7; , a card for the coefficient of 
that term is included in the tabulation ahead of the summary cards. This 
automatically adds the new coefficient to each term of the series. When the 
subscript 7 in 77’; changes, the new coefficient card must form a separate total; 


1 This operation is generally known as progressive totalling in machine operation. 
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when it does not change, the coefficient card must tabulate in the first summary 
card total. 


To illustrate the tabulation of power tables, the formula for the cube table is— 
3 14 2rn 
¢ =6 724+ °7T). 
The successive operations yield the following table: 
TABLE III 


Operation number 


3 


0 1 
6 7 
12 19 
18 37 
24 61 125 
30 91 216 
36 127 343 
42 169 512 
48 217 729 
54 271 1000 


oe oOonNnoaoartk WN 
DADAAAAAAAS 


— 





In actual machine work, operation 1 can be omitted and work begun with opera- 
tion 2. The machine is set to add the coefficient 6 of the highest term from 
each card and an accumulated total is printed and punched for each card tabu- 
lated, giving the results shown under operation 2. An additional card is punched 
for the coefficient of the second term, 1, and placed ahead of the cards produced 


in operation 2. The cumulation and punching is repeated, giving the results 
shown under operation 3. The summary cards from this operation are cumu- 
latively tabulated, giving the results shown under operation 4, which is the 
table of cubes desired. 

Similarly, for a table of the fourth power, the formula x* = 24°7; + 12 ‘7, + 
2°T, + °7; indicates the following operations— 


TABLE IV 


Operation number 


3 








0 
14 
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Note in operation 3 where the subscript does not change, the coefficient 2 is 
added to the first card punched by the machine, while in operation 4 where it 
changes, the coefficient 1 appears as a separate total. 








4. Tables of polynomials. To tabulate values of f(x) = a + br + cz’... 
(where a, b, c, --- , are positive or negative coefficients) the method is similar 
to that of preparing power tables except that the coefficients to be added are 
determined by multiplying the coefficients of the formulas for the different powers 
by the values a, b, c etc., adding the coefficients of like terms in the various 
formulas, and using these resultant coefficients in place of the simple coefficients 
used in the power tables. Thus if we wish to tabulate values of f(z) = 4 + 32 + 
2x” + x’ the coefficients are found as follows: 













4o® = 4'T, 
+32 = +377, 
+ 22° = + 2°T, + 2.2°7, 
+ 2°= + *T, + 30‘*T, + 120 °7; 





fa) = 44 OT + ATs + 30°Ts + 120°T, 


This equation gives instructions to perform six operations with 120 as coeffi- 
cient; adding the coefficient 30 as a separate total when there are 4 operations 
remaining; adding 4 to the first summary card total when there are 3 operations; 
adding 6 as a separate total when there are 2 operations remaining; and adding 4 
on the last operation. 

The first few totals appear thus— 


TABLE V 








Operation number 
4 
0 4 
1 0 0 0 0 6 10 
2 | 0 0 30 34 40 50 
3 120 120 150 184 224 274 
4 | 120 240 390 574 798 1072 
5 | 120 360 | 750 1324 2122 3194 
6 120 480 1230 2554 4676 7870 
7 | 120 600 1830 4384 9060 16930 
8 120 720 2550 6934 15994 32924 
9 120 840 3390 10324 26318 59242 
10 120 960 4350 14674 40992 100234 

















It is not necessary to confine these tables to values for whole numbers, as we 
can tabulate equally well values of f(x) for intervals of x of .1, .01 or .001 or 3, 
3, z etc. In this case, before combining formulas for different powers we multi- 
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ply both sides by the desired interval raised to the power to which z is raised in 
that particular formula, then add like terms as before. 


To tabulate the previous example in .1z intervals we proceed as follows: 
4x° = 4.000 ‘Ty 
32/10 + 3°T 
2(2/10)" = + .02°T, + .04°7; 
(x/10)° = + .00001 *7, + .00030 ‘7, + .00120 °7; 
f(z) = 4'°To + .32001 °7, + .04°T, + .00030 ‘7, + .00120 °7; 
TABLE VI 


Operation number 
- 6:f(z) 


0 4.32001 

0403 4.68032 

0418 5.08243 

0457 5.53024 
| 0532 6.03125 
| 0655 6.59776 
| .0738 7.23807 
| 0993 8.07768 
| 1332 8.95049 
| .0096 =—.0435 1767 1.04951 1000000 


1 
2 
3 
4 
5 
6 
7 
8 
9 
10 





Where any coefficients are negative in the equations expressed in ’7; terms, 
they are simply added in as minus figures. 

To round off the preceding function to 3 decimal places, we add 5 to the con- 
stant term '7') in the position to the right of the last decimal retained, i.e. in 
this case the 4th decimal place. The constant term is then 4.0005. 


Exact Counter reads Prints 
4.32001 4.32051 4.320 
4.68032 4.68082 4.680. 
5.08243 5.08293 5.082 
5.53024 5.53074 5.530 
6.03125 6.03175 6.031 
6.59776 6.59826 6.598 
7.23807 7.23857 7.238 
8.07768 8.07818 8.078 
8.95049 + 8.95099 8.950 

10.00000 ; 10.00050 10.000 


5. Automatic calculation of polynomial coefficients. Frequently when 
polynomials are being evaluated, the process of forming the coefficients can be 
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performed automatically from a punched-card table. Such a table consists of 
set of cards for each power x’ containing the multiples of all the coefficients of 
each of the terms ’7’; in the formula (3) for that power. These multiples are 
1, 2,3,4, --- , 9; 10, 20, 30, 40, --- , 90; 100, 200, --- , 900; 1000, 2000 etc., 
and may be produced automatically by making a linear table of each coefficient. 
in the manner described in this paper. Each card is punched with the informa- 
tion called for by the heading of the following card form: 


4 | 
/ 8 j i multiple 


coeff. x 
multiple 
07 06 03 00005 008400 


The particular figures indicated are those which would be punched for the 
term 5(1680)°7; in the representation of 52’ according to formula (3). 

The table is used by withdrawing the cards for the coefficients a, b, c, d, ete. 
’ of the desired polynomial. For instance, if one of the polynomial coefficients is 
14485 2’, we select from the z’ section of the table all cards containing the multi- 
ples 10000, 4000, 400, 80, and 5. In the z’ table there are 4 cards for each multi- 
ple, one each for terms °7,,°T;, “72, and’7,. These cards are combined with 
the cards selected for the other coefficients of the polynomial and sorted to bring 
all cards for each ’7'; together. The cards for each term ’7'; are then automati- 
cally added on the electric accounting machine. 


6. Subdividing tables. In preparing tables it may be desired to prepare 
the table in more detail at certain points, giving values of the function at 1/10, 
1/20, 1/50, or 1/100, etc., of the interval of the rest of the table. This may 
readily be done by recalculating the coefficients of the cumulative terms, and 
using these values in the same manner as the original ones. 

There are many formulas for the determination of the subdivided differences 
given in various texts on interpolation, such as those given by Comrie [8] and 
Bower [9]. One effective method is to use formulas (3) to calculate the sub- 
divided differences. The values called for in the formula for the highest power 
are taken from the table of the function at the regular interval, giving effect to 
the rule involving subscripts. These coefficients are reduced by an amount 
sufficient to cancel the coefficiént of the highest cumulative term, and the coeff- 
cients of the remaining cumulative terms are reduced in proportion according 
to formula (3) for the highest power. Usually the coefficient of the highest term 
of the formula will divide evenly into the coefficient taken from the table, and 
the other reductions are calculated by multiplying this result by the other 
coefficients of the formula. The highest remaining coefficient is then reduced 
by an amount sufficient to cancel itself, and, by use of the formula (3) for the 
power whose highest cumulative term matches the highest remaining coefficient, 
the reduction to the remaining cumulative terms is calculated and subtracted. 
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The highest remaining coefficient is reduced in a like manner, and this process is 
continued until all the cumulative coefficients have been analyzed. 

The partial cumulative coefficients thus computed are multiplied by the de- 
sired subdivision 1/m raised to the power of the corresponding formula (3), 
and recombined to form the new coefficients, as shown in the example below. 
In taking values from the table, when the subscript does not change, the tabular 
value must be reduced by the amount of the higher coefficient with the same 
subscript, to give effect to the rule that the coefficients in such cases are incre- 
ments (see last example in section 3). 

To subdivide the polynomial of section 4 at z = 7.0, we take the italicized 
values from Table V starting at f(7) as 'T) , and proceed as follows: 











6 Te “7 4 Te 3 Ts 2 7, 1 T, 
960 10324 
From Table V ..... 120 -120 3390 -3390 15994 16930 





840 6934 


















840 6934 
840 420 70 35 


6864 









6864 
6864 


If the interval is 1/10 we have: 





6 T's . i "Ty. 2 7, 1 T. 
2°/10° = ,00120 .00030 .00001 
352/108 = 0840 .04200 .0070 .00350 
4902°/10° = 2.94000 .49000 
34322°/10° = 68.6400  34.32000 
120362/10 = 1203.60000 + 16930 


f(x) = .00120°T; + .0840 °7; + 2.9823 ‘7, + 68.6470 *T, + 1238.41351 *7, + 
16930'7) provides the coefficients for subtabulating the function at the desired 
interval, beginning at the argument x = 7.0. 





1. Accuracy of Tables. When the values of the coefficients are not 
exact, owing to the original values for a, b, c etc. or the dropping of decimals in 
the computation of the coefficients, the errors accumulate fairly rapidly. Each 
coefficient will introduce its own error into the summation. 
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To maintain accuracy throughout a long table it is advisable to transform f(z) 
by Horner’s method of decreasing the roots [10, pp. 100-101], compute _new 
coefficients for the transformed equation at intervals, and prepare the table in 
sections. Decreasing the roots by r gives us a new starting point at x = r. 

Since two or more functions may be computed at one time, a. function for 
which the coefficients are not exact may be computed by adding in the usual 
way from the starting values and subtracting from the ending values simul- 
taneously. As many digits as agree in both tabulations of the function may be 
considered correct. 

The tabulations can be made to practically any degree of accuracy on the 
equipment available, as the newer machines can be formed into counters of any 
capacity up to 80 digits. In practice, counters of 16, 20 or 24 digits will ordi- 
narily suffice for the accuracy desired and two or more functions can be evaluated 
simultaneously. Cards are read and added at the rate of 150 per minute, or read, 
added and listed on the tape at the rate of 80 per minute and new summary 
cards produced at the rate of 40 per minute (on alphabetic equipment with gang 
summary punches). Computation may be carried out with additional decimal 
places and the final tabulation of the function rounded off to the nearest number 
retained. 


8. Summary. “The cumulative or progressive-total method is shown to be 
applicable to the preparation of tables of functions expressed in the form of 


@ power series. 

The cumulative formulas for the powers through the twelfth power have been 
presented, and simple methods are given for transforming a power series into its 
corresponding cumulative formula, for changing the interval of the table, 
rounding off the values of the function, and subdividing the table at desired 
points. 

It is hoped that this discussion will make tables in printed or punched-card 
form more generally available as a tool for the computer. Since tables may be so 
readily prepared by this process, the usefulness of the tabular method of solving 
problems is greatly increased. 

The author wishes to acknowledge his thanks to Professor P. S. Dwyer for 
various suggestions, particularly in connection with section 2. 
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Introduction. Let F,,---,E,, denote mn _ arbitrary events. Let 
Doi ...iv¢4,---9; , WhereO S$ 7 Sj S nand (y, ---, v;) is a combination of the 
integers (1, --- ,n), denote the probability of the non-occurrence of E,, , --- , E,, 
and the occurrence of £,,,,,---, E,;. Let py,....,; denote the probability of 
the occurrence of E,, , --- , H,, and no others among the n events. Let S; = 
Zp,,....; Where the summation extends to all combinations of 7 of the n integers 
(1,---,n). Let pa(m,---,), (1 S m Sk S n), denote the probability of 
the occurrence of at least m events among the k events E,,, --- , E,,. 

By the set (1, --- , %, +++ ,2a) — (41, ---, 2%) (where b S a) we mean the 
set (241, °-+, 2%). And by a E -combination out of (7, --- , Za) we mean 
a combination of b integers out of the a integers (11, --- , Za). 

We often use summation signs with their meaning understood, thus for a fixed 
k, 1S k S n, the summations in =p,,...,,, or 2pm(vi, --- , v%), extend to all 
the (7 }combinations out of (1, ---, n). 


The following conventions concerning the binomial coefficients are made: 


0 a ‘ , 
(3) = 1, (5) =o if a<b or if b6< 0. 


It is a fundamental theorem in the theory of probability that, if Z,,---, En 
are incompatible (or “mutually exclusive’’), then 


Pi(l, +++, ) =~Mt eee + Da. 


When the events are arbitrary, we have Boole’s inequality 


p~(l, +++, 7) Spat rare + Dn. 


Gumbel’ has generalized this inequality to the following: 





) < Zpi(r1, Pine » Vk) 


cp 


rll, ween 


1C. R. Acad. Sc. Vol. 205(1937), p. 774. 
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fork = 1,---,m. The case k = 1 gives Boole’s inequality. Fréchet” has 
announced that Gumbel’s result can be sharpened to the following 


(1) age = SBM ==> 2s) 5 Beis --- 1 
n—1 n-—1l1 
k k-1 
fork = 1,---,n-— 1. Thus, A; is non-increasing for k increasing. On the 
other hand, Poincaré has obtained the following formula which expresses 
pi(l, «++, %) in terms of the S,’s, 
pl, ow’ n) — Zz Pr, — Loin + 2. Posters = ©e 


+ (—1)"preen = . (-1)*"'g,. 


A, 


(2) 


In the present paper we shall study the more general function pm(v, --+ , ve) 
as defined above. First we generalize Poincaré’s formula and Fréchet’s inequali- 
ties. In Theorem 1 we establish (for 1 < m S n) 


Pally +++ 50) = DPreine — (7) Dparernes 


+( : ') Do Pressman Feet (—y(" 


— (-»'(” ahi " ia, 


t=0 
Although this result is well known, we prove it in preparation for Theorem 2, 
Theorem 3 establishes 


(4) Am = ZPm(Vi, +++, Meta) © VPm(vi, +++, Me) _ A” 
? 


k+l = = — 
n—m n—m 
fds ual _— 


fork = 1,---,n—1landl Sm S.-k. 
Next, we extend the inequalities (4), and in Theorem 4 we show that 


(5) Ay” S 3(Ai) + AR) ; 


which states that the differences A; — Aisi (Kk = 1, --+ , m — 1) are non-decreas- 
ing for increasing k. From this and a simple result we can deduce (4). Also 
Theorem 2 establishes that 


@) Da FF) Ses pall, ms (371) San, 


t=O 


? Loc. cit., Vol. 208(1939), p. 1703. 
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for 2] + 1 S n — mand 21 < n — m respectively. These inequalities throw 
light on formula (3) and are sharper than the following analogue of Boole’s 
inequality for p,(1, --- , »), which is a special case of (4): 

(7) Pm(1, +++ 2) SF LPy...m - 


The last statement will be evident in the proof. 
In Theorem 5 we give an “inversion’”’ of the formula (3), i.e. we express 7)... 
in terms of the pm(v , --- , v%)’s, as follows: 


- 1 
a = D> Pm(v1, a — > Pm(r1, 22+) Vm+i) > 20% 


(8) + (-1)" pm(1, +++, 2) 
7 (—1)' 2. Dm(r1, Sees Vm-+i)- 


t=0 


This of course implies the following more general formula for pa,..-a, ; 


—]1 r—m ‘ 
. ss ) Pa;:++a, = a (—1) >. P(r, oes Vm+i) 


where (ai, --+ , a) is a combination of the integers (1, --- , n) and where the 
second summation extends to all the (,, . ; combinations of (a1, +++, a). 


Since it is known’ that we can express other functions such as S, , i eee Oe 
terms of the p,,....,’8, we can also express them in terms of the pm(v, - ++ , vx)’s, 
provided r = m. 


Finally, for the case m = 1, we give in Theorem 6 an explicit formula for 
Pu...r) in terms of the pi(v1, --+ , v%)’s, as Shown in (9), 


Puen = — pir +1, +++, n) + DV piln,r $1, «++ , 2) 


— Di pili,ve,r+i,---,n) + oes 


¥1»¥Q 


+ (-—1)" > pil, :-- »r,r+l1,---,n), 
= 2 (-1)™ » pir, Se ee i 1, oN), 


where (»,,---, »:) runs through all the (7)-combinations from (1,---, 7). 


This of course implies the following more general formula: 


- 
Play-+ ay] = -™ (—1)*" 2. pilri, ee » Yi, Ar+1, 2 , Qn), 
i=l (v1.5 * og) 


+4 Fréchet, “Condition d’existence de systemes d’événements associés 4 certaines 
probabilités,’’ Jour. de Math., (1940), p. 51-62. 
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where (a1, °-+,Q-,-+:+ Qn) iS a permutation of (1,---,n) and where 
(vy, +++ , ¥) runs through all the (*))-combinations out of (a,,---,a,). From 


Theorem 6 and two lemmas we deduce a condition of existence of systems of 
events associated with the probabilities p:(y,, ---,vm). The author has not 
been able to obtain similar elegant results for the general m. Probably they 
do not exist. 


2. Generalization of Poincaré’s formula; Generalization and sharpening of 
Boole’s inequality. 
THEOREM 1: 


a, (") ae 
+ (” ? ') DL Pr--tase + (ayn (" 7 ) Pr 


n= 


(3) 


Proor: We have 
(10) Pn(1, i, n) = d Zz Plur'**um+ol > 


where the second summation extends, for a fixed b, to all the ( ) combina 


n 
m+b 
tions of (1,---,%). Further we have 


n—m—C 


(11) Pye °— > stint 


where the second summation extends, for a fixed d, to all the an - elie 
combinations of (1, --- ,n) — (4, -++, mse). The formulas (10) and (11) are 
evident by observing that the probabilities in the summations are all additive. 
Now we count the number of times a fixed pjy,...u,4,) appears in (3). -By (11) 
this is equal to the sum 


Ye _ (m\(m +b 4 m+1\(/m+ b\ _ 
m 1/\m+1 2 m+ 2 _ 


Herta ina) 


since this number is the coefficient of (—1)”x™ in the expansion of 


(—2)"* (1 i sy" = (-1)"2"(1 — 2)’. 


Thus by (10) we have (3). 
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THEOREM 2: For 21 S n — mand 21 S n — m respectively, we have 


2+ 1\ au - 
6) XY (- ot Y) Susi S Pall «= n) So (—1yi(™F YBa 
i=0 u F i=0 t 


Proor: By the reasoning in the previous proof, it is sufficient (in fact also 
necessary) to show that 


f/m —1+i\/m+ob f/m —1+i1\f/m+b 
> ( a ce a | a mae) <a, 


Since 
e —1+ _ + ‘) _ (m+)! (*) 1 
i m+i) (m—1)!Db! m+i 


is an integer, it is sufficient to show that 


2U ‘ h 1 21+1 ; b 





IA 
© 


m+2 


Suppose b > Oiseven. For7z Ss b/2 — 1, we have Ti > 1 so that 
7+2 mt+i  st+l1 


Also form = 1. Hence 


{41° mt+it+l1l~ i+2 - 
( b ) 1 _b-t m+i (*) 1 
itij/mt+it¢1 titimt+i+ti\i/m+i 


eH Qan Oa 
~i+1i4+2 m+i \VJm4+i 


; | 1 so that . ae Ft. < land 


z+im+i7+1 


( b ) 1 5 <(j) 1 
ti+1l/m+i+1 tiJ/m+ia 


Thus the absolute va'ies of the terms of the alternating series 


b! 
y(- in) (m + b)!(m — 1)! 


are monotone increasing as long asi S —~ — 1, reaching maximum at7z = 3 and 


« 





2.) oO 
~ 
IV 


Ti | 




















For i = b/2 h : 
ori = b/2 we -" iz 





nl o 


then become monotone decreasing. 
Therefore (12) evidently holds for 21 < b/2 and 21 + 1 < b/2 respectively. 


For t = 






: + 1 we write 


b! - i(b 
y(- (oh (m + b)!(m — 1)! ~ 2 1) () an m+1 


' b! _ b—t—1 - j b _ -. 
“etumcm- & () anes 








nv 


id 
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b! 1 
a 
m+ bdim—D! = m+ , We see that the 
righthand side is an alternating series whose terms are non-decreasing in absolute 
values. Hence (12) is true. 


If b is odd, the case is similar. 


From the above and the fact that 


3. Generalization of Fréchet’s inequalities and related inequalities. Before 
proving our remaining theorems, we shall give a more detailed account of 
the general method which will be used. In the foregoing work we have al- 
ready given two different expressions for the function p,,(1, --- ,), namely, 
formulas (3) and (10), but they are not convenient for our later purposes. 
Formula (3) is inconvenient because it is not additive and because the p,,...,,’8 
are related in magnitudes; while formula (10) has gone so far in the separation 
of the additive constituents that its application raises algebraical difficulties. 
Let us therefore take an intermediate course. 


Let each (”)-combination (v1, -++, Mm) out of (1, --- , n) be written so that 


vy, < ve <+++< vm. Then we arrange them in an ordered sequence in the 
following way: the combination (v1, ---, vm) is to precede the combination 
(ur, °° , Mm) if, for the first », ~ u;, we have vy; > w;. After such an arrange- 
ment we symbolically denote these combinations by 


was [(*)]. 


Further, all the (*) combinations out of (4, --- , %) where the latter is a com- 
bination out of (1, --- ,m) are arranged in the order in which they appear in 
the sequence just written. For example, all the G 
(1, 2, 3, 4) are ordered thus: 

(12) (13) (14) (28) (24) (84). 


) combinations out of 


Let U denote a typical combination (4, ---, um). By Evy we mean the com- 
bination of events E,,,---, Eu, so that pu = p,,...,,- In general, let the 
combinations U,, --- , Urs, U's be given, then pu;...uj_,v, denotes the proba- 
bility of the non-occurrence of U7, , --- , Us+ and the occurrence of U,. 
Now let I, I, ---, | (*) _ | = J, [@) = Z denote all the (*) com. 
m m m 
binations out of (vy; , --- , ») in their assigned order. We have 
(13) Poli, +++ Ve) = Prt Pru + Privm + +++ + Pr..vz. 


This fundamental formula is evident. Of course it is possible to identify the 
p’s on the right-hand side with the ordinary p,;...,;’s, but we shall refrain from 
so doing and be content with the following example: 


po(1, 2, 3,4) = pro + piers + Prorsa + Pres + Prasa + Prom - 
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THEOREM 3. Fork = 1,---,n—1land1 Sm Sk we have 


E = ™) lpm, ee s (, se Lpm(r1, 22+, De). 


Proor. Substitute (13) and a similar formula for k + 1 into the two sides 
respectively. After this substitution we observe that the number of terms is 
the same on both sides, since 


(em) + OR) = (ea 7 En): 


Also, the number of terms with a given U = (yu, ---, wm) unaccented is the 


same, since 
n—m n—-m _ n—m n—m 
k—-m/J\k+1—m)/ \k+1—m/\k—m/)’ 


Let the sum of all the terms with U unaccented in the two summations be 


denoted by ox41 = on%41 (ur, +--+ , Mm) aNd o, = ox (41, --- , Mm) respectively. It 
is sufficient to prove that 


as (t= mows (et ir) 


n 


for any U. o; contains c “a _ terms each of the form 7,;...,;4,-..u, Where 


OslSun—m and where (1 oe ee Mm) is a Lm })-combination 


out of (1,---, wm). For fixed (u:,---, um) and a fixed J but varying 2’s, ox 

contains k = a 1 terms of the form 7,;...1iy,---um » With exactly | accented 

subscripts. Let the sum of all such terms be denoted by of”. Evidently of” 

nn (*°  ™ 
l 


) terms. As a check we have 
n — Bm Um — ™ nr — Um Km — ™ 
Geo GT) 
2 — tm\({tm —m\ _ (n—m 
+ - )=(¢=%) 


which is the total number of terms in o; . 
We decompose these p’s partially, as follows: 


Bm—m—l 


Poj-+-¥imi-+ +m Pri -+-¥itcmi:*-Bm+d? 


b= m+ie- + om +d 
where (1, +--+, Vine, Mi,» *** » Mm4d) IS @ permutation of (1, --- , um) and where 
the second summation extends, for a fixed b, to all the + es . ~ *) combina 


tions out of (1, --+ , um) — (¥1, +++, mt, Mi, °°? Mm). 
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Now consider a given 
Doj--:ptdyee-Aghi s+ *m 
where 0 S t S um — m and (1 -->+ parr -++ Agi +++ Mm) iS a permutation of 


(1,+++, um). It appears (;) times in of”. Hence it appears 


(Go) tert G)+ GRE = Cae) 


times in o . 
Therefore to prove (14) it is sufficient to prove that 


Ce < ws. 
lk -—m k+1—m/>\k+1—m™ k—m ‘ 
By an easy reduction we have 


(n—pmti—-k+m)sn-k 


— pa tt+msd; 


since t S um — m this is obvious. 
THEOREM 4: For2 Sk ZSn-—1land1 sm & k we have 


LPm(r1, soe, Me) 1 Lpm(v1, “= » Ve-1) 1 Lpm(r1, oy Mest) 
pS og 5 ee et 5 5 
n—m 2 n—-m 2 n—-™m 
k—m k-—-1l—m™ k+1—m™ 


Proor: By the reasoning in the previous proof, it is sufficient to show that 


9 oa n—m oe 
k-—-1—m/\k+1—m™ k—m 
n—m n—m n— ta tt 
=\k —m/\kK+1—m/\k -—1l1—m™ 


4{?-* n—m™ Nn— pn tt 
k—-m/\k -—1—m/\k+1—m/’ 


for0 S$ t S um — m. By an easy reduction this is equivalent: to 
2(n — k)(n — pm tt —kK+m+1) S (n—k+4+1)(n — k) 
+ (n— pm tt—k+m+1)(n — un tt—k+m) 


(n — um +t —k+m+ 1)(um — t — m) S (n — k)(um — t — m). 
For t = um — m we have equality, otherwise we have 


—tm ti+tm+1s0. 
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We can deduce Theorem 3 from Theorem 4 and the following result (a case 
of generalized Gumbel inequalities): 










. — 1 
(15) (” ae >) Pall Po Te n) = lpm, ee » Med 


Proor oF (15): Substitute from (13). Consider the p’s with U unaccented. 
The number of such terms is the same on both sides. But on the left-hand side 
they are all the same pr’...(v—1’v , While those on the right-hand side, being of 
the form py;...vju Where 0 S AX S U — 1 and (U,, --- , Uy) is a combination 
out of (1, --- , U — 1), are greater than or equal to it. Hence the result. 


4. The puz,.--a,’S in terms of the p,,(v, - 


of the pi(vi, ---, %)’S. 
THEOREM 5: Forl S m S n we have 


—1 
(” seal 1) Pa De Pm(P1, oe Vm) mod De Pm(r1, Se Vn+1) + .:-- 


-, v)’s and the pja,---a,)’S in terms 








(8) + (-1)" "pn(1, ---, ) 
= 2 (—1)' DD p(t, +++ Hts) 


Proor: As in the proof of Theorem 3, consider ox(ui,---, um). Here 
msxk«zn. Since a given 


(16) 






Poi-++ptda-+-Agit * "Hm ’ 





appears s~ tn + 
Pp} k—m 


eave. if n—pmwt+t21, 
J :. if n— um +t = 0. 
times on the right hand side of (8). Hence for fixed (u,--- , um), the only 
p’s of the form (16) which actually appears are those with t = um — n. But 
lim Sn, thust = 0, um = n, and (Ay, --- As, Mi, °** 5 Mm) iS & permutation of 
(1,---, mn). The term in question is therefore p,..... Since the number of 
("" ) combinations of (1,---, n) withun = nis (” ei 7 , we have the theorem. 


THEOREM 6: For 1 Sr S n — 1, we have 


) times in o,, it appears 








Pu---r]) = — pilr + 1, rr») +> piln,r+1,---,n) 
(9) — DL pli,me,rti,--+,n) +--+» +(-1)" D pill, --- , m) 


~ 2, (27 7 pin, *** mF +1, +++, md, 
t= P40 ME 
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where (v1, +++, vi) runs through all. the ({)combinations out of (1,---, 7). 


Proor: We rewrite (14) for the special case m = 1, 


(17) Pili , ie » Mk) = Pu + Pui ue i Pui ---wi-iee > 


where wi < we <-+- < yx. Substitute into the right hand side of (9). After 
the substitution let the sum of all those p’s with uw unaccented be denoted by 
o,. The terms in o, are of the form p,;.../_,. Where 1 Ss S uw and 
(ui, *** » Ms—1) is a combination out of (1, --- , uw — 1). 

First consider a fixed » S r. For a fixed p,;...,:_,, we count the number of 


times it appears in o, , that is, on the right hand side of (9). This is evidently 
equal to 


: oo r—pte Yo 0, ifr—p2il, 
“ea 
2 ( j~&8 2X | : j—é& 1, ifr—p=0. 


Thus the only terms that actually appear are those with u = r; and each of such 
terms P,'...u;_,, appears exactly once with the sign (—1)°. Hence their total 
contribution is 


(18) Pr — zz Poir + } Pa tos + (—1)"" py... = Pi---r, 
¥i+¥2 


v1 


by an easy modification of Poincaré’s formula. 

Next consider a fixed » 2 r+ 1. Every term with » unaccented in ao, is of 
the form (with the usual convention for uw = r + 1) Dyl...us 4ny’-.-w—’s » Where 
(u1, -*+, Ms) iS a combination out of (1, --- , 7); and it appears exactly once 
with the sign (—1)*. Their total contribution is therefore 


— Detn’---@—'a + z, Poy (r+l) "+++ (u—1) pp 7 Pojivg(r+l)’ +++ (u—l) p eee 
v1 


¥1¥2 
r—l 
ced (—1) Pi'++- (ut) "pp = — Pi--er(r4))’+++Gu—l) "py 


by another application of Poincaré’s formula. Summing up for yu = 
r+1,---,n, we obtain 


(19) — (Pr...r(r41) HF Di...r(rgty (rg) Hees + PA. -r(r-41)? ++-(n—2)'n) 
Adding (18) and (19), we obtain as the sum of the right-hand side of (9) 
Pi... — (Pr..-r(r4-2) + Dh...r(rgiy(rga) Hees + Di..r(r$1)"--(n—2)!n) 


= P1..-r(r+1)! (r42)! en? = Pii---r] 


by an easy modification of (17). 


5. A condition for existence of systems of events associated with the proba- 
bilities pil pire y Vy). 
Lemma 1: Let any 2” — 1 quantities q(a:,---, ax) be given, where k = 
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1,-++, , and for a fixed k, (a1, +++ , ax) runs through all the : -combinations 
out of (1,--+,n). Let the quantities Q(a , --- , ax) be formed as follows: 
Q(0) = 1 — q(l, ---, n), 
Q(a1, ree, a) =— q(ars1, eee , An) + >» an, Oks1, °°"; an) 
_— a gin, Vo, Aki,» °°", Qn) oe coe + (—1)** (1, vee n), 


Vie 


where (v; , --- , v:) runs through all the (‘)-combinations out of (1,---,n) — 


(Qni1,°°*, Qn). Then the sum of all these Q’s is equal to 1. 
Proor: Add all these Q’s and count the number of times a fixed g(u1, «++ , ux) 
appears in the sum. For 1 S k S n this number is equal to 


=" (‘) * (*) end —y(7) = 0. 


Hence we have the lemma. 

LremMaA 2: (Fréchet) Given 2” quantities Qta,...a,) Where (a1, +++, @) runs 
through all combinations out of (1, --- , n) including the empty one. The necessary 
and sufficient condition that there exist systems of events E,, --- , E, for which 


P{a---a,]. = Dont 
(where po; denotes the probability for the non-occurrence of F,, --- , E,) ts 
that each Q = 0 and that their sum is equal to 1. 
Proor: Since the probabilities pja,...c,; are independent, i.e., unrelated in 
magnitudes except that their sum is equal to 1, the lemma is evident. 
THEOREM 7: Given 2” — 1 quantities q(ai, --- , ax) as in Lemma 1, the neces- 


sary and sufficient condition that there exist systems of events E, , --- , En for which 
pila , sas » Ak) a g(a , — » Qk) 
is that for any combination (ar41, +++ ,@n),1 Sr Sn — 1, out of (1, --- , n) we 


have 7 


— glaris, 7° Gn) + Dy qa, » Ors; °°) an) v= i Ql ar, 5 Ong » Or4ty °°, aa) 
v1 ; 


Vive 
+ an + (—1)""* (1, 7295 n) = 0, 
and thus 
_- q(, Pry 


Proor: The condition is necessary by Theorem 6. It is sufficient by Lemma 
1, 2 and an obvious formula expressing pi(a; , --- , a) in terms of the py,...»,)’s. 


,n) 20. 
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NOTES 


This section is devoted to brief research and expository articles, notes on methodology 
and other short items. 
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As far as the author is aware, H. C. Carver was the first to point out that 
while the formulae ordinarily given for Sheppard’s corrections for central mo- 
ments are valid for moments computed about the population mean, there are 
still systematic errors present when they are applied to central moments calcu- 
lated from any particular grouped frequency distribution [1]. This is due, of 
course to the fact that the mean of a grouped frequency distribution is in general 
different from that of the distribution before grouping. For a fixed class interval 
k, Sheppard’s corrections give the average value of a moment about a fixed 
point of a given order for all the groupings of this class width possible and will 
fail to do so if the moment in question is calculated for each position of the class 
limits about a point which varies as the class limits shift. Thus Carver [1] 
pointed that the commonly used formula (for a continuous variate), 






k? 


(1) fom 8 ™ Ta 





should, if v. is calculated about the mean of the grouped distribution as it is in 
practice, be replaced by 


k? 
(2) nano ton 


in which oy is the variance of the means of grouped distributions over all posi- 
tions of the class limits with the fixed class width k. 

Recently J. A. Pierce [2] gave a method for deriving the required formulae of 
the type of (2) and gave actual formulae for both moments and seminvariants 
through the sixth order. It is the purpose of this note to point out that the use 
of moment generating functions provides a more elegant and concise way of 
arriving at formulae equivalent to Pierce’s though in a somewhat different form. 
This method can be immediately extended to distributions of two or more 
variates. 

In a previous paper [3] on Sheppard’s corrections for a discrete variate, the 
author made use of the following argument: It is assumed that for a fixed class 
width k, any point in the scale on which the variate z is plotted is as likely to be 
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chosen as a class limit as any other; choosing a system of class limits for grouping 
the data is then equivalent to placing at random on the z-axis a scale with 
division points at intervals of k. Once the system of class limits is chosen any 
value of xz before grouping bears to the class mark, 2x; , of the class in which it 
falls the relation, 


(3) V=xrte, 













in which z and e¢ are independent variates. The frequency law governing 2, is, 
of course, that of the population from which it is drawn while « is distributed 


in a rectangular distribution with the range (- + :) for a continuous variate 
m—1,m™ 1 . ' ‘ , 
and | — —.—— k, ——— k) if m consecutive values of a discrete variate are 
2m 2m 


grouped in each class interval. In either case 
(4) 


in which M,,(#) is the moment generating function of the variate z;, etc. The 
expansion of both sides of (4) in powers of # gives the relations between the 
average values of moments of the grouped distribution over all positions of the 
scale and the moments of the ungrouped distribution from which Sheppard’s 
corrections are obtained by solving for the moments of the ungrouped distribu- 
tion. The relations are valid for any fixed point about which the moments are 
computed; if this fixed point be taken as the mean of the ungrouped distribution 
the ordinary Sheppard’s corrections for central moments result. 

But it is quite easy to modify (4) to give the necessary relations in case the 
moments of each grouped distribution are computed about the mean of that 
distribution. We have only to write 





M;, (#) = M,(3)M (#8) 













(5) GZ=2%-EF+z 
in which Z is the mean of the grouped distribution for which 2; is one of the class 
marks. Then 





M.,(3) = M,,-2,2(0, w) Lent 











(6) 


If we write, 


Arask;—2,2 = Are ? 





in which X,, is the product seminvariant of order rs of moments about the means 
of the grouped distributions and of such means, the expansion of the logarithm 
of the second member of (6) gives 





- (X10 + Xo)” 3 : oe, 





(7) 1 + (Aso + Xor)d + (Avo + 2An + Xn) 
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in which 


(Aro + An)” — dw + Try-1,1 + ++. + 4 Ark yk + eee + Nor . 


The expression of the logarithm of the right member is’: 


3 3° 2 4 B,k** 1 3” 
(8) MO + Mat Mate + dy (8 ( _ Js) Gav 
for a discrete variate (the result for a continuous variable is obtained merely by 
letting m — ©) in which }, is the rth seminvariant of the ungrouped distribution 
and B, is the sth Bernoulli number. 
We may without loss of generality take the origin for z at the mean of the 
ungrouped distribution so that 4, = 0. Further it is easy to see that 


\ir = 0, r=0,1,2,3,---. 
Consider 
For a fixed Z, i.e., for a given grouping, this becomes 
E(x; — 2) = 0 
Then since »;, is the average of this over all groupings with a given class interval, 
i, = 0, and from the expression for \i, in terms of the moments 3;; it is obvious 
that also A1, = 0. 


Then we must also have Xo = 0 as is otherwise obvious and (7) can be rewritten 
bs _ . ¥ " . _. 3 
(9) 1 + (ao + oe) = + (Aso + BAar + dos) 55 + ++ 


Now from (8) and (9) by equating coefficients of like powers of 3, we get the 
set of formulae: 


A = 0 
. i 1\k 
Ae iw + Bn — (1 - 2) 
Aso + Bre + Avs 


hee + din + Oi + ha + (1 - 2 


These formulae, however, do not give the sought Sheppard’s corrections for 
seminvariants calculated from grouped distributions of a discrete variate. See 
below. 

Referring to formula (10), p. 58 of the author’s paper cited [3], it is easily seen 


by comparison that the required moment formulae are obtained from the general 
formula 


[n/2] 
(11) ie = > (;.) 024(F10 + we, 
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in which ae, is given by formula (9) of this former paper. For = 1, 2, 3, 4 
we write down immediately 
wm = 0 (D1 = 1 = O) 
lal ( 1 
Me = do + io — 1--,; 
m? 


v3 + 3¥e1 + Hos 
= v9 + 4931 + Gi20 + Dos 


i. ole we 3\ k 
age (1 _ 1) (d20 + Hea) 9 i (1 LY) 3) x 


In these formulae, 7,9 is, of course, the average value of rth central moments 
about the means of grouped distributions. From the definition 7,.(s # 0) is the 
average value of the product of the rth central moment of a grouped distribution 
by the sth power of the mean of the same grouped distribution. Also, it must 
be noted that in the formulae (10) the X,.’s there are to be calculated by the 
usual formulae from the moments, 7;; , and are not themselves the average values 
of like seminvariants calculated from the separate grouped distributions. Thus 
though the formulae (12) give the sought Sheppard’s corrections for moments, 
the formulae (10) do not do the like for seminvariants in general. However, 
since in each grouped distribution, 


de 
and 
As = Bs 
we have, taking the expectation or average value over the grouped distributions, 
E(d2) = E(v2) = 29 = Yeo 
and 


E(d3) = E(vs) = t90 = Avo, 


and the first two formulae of (10) do give the Sheppard’s corrections for A: and 
3 calculated from grouped distributions of a discrete variate. 
But the case for , is different. In each grouped distribution, 
N = %4-— 373 , 


and if we define l, by 
E(\,) = 4, 
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we have 
m= 3E (v3) 
= i — 3 (P20 + V2:y5) = ao = SV 2:05 ’ 


if v;», is the variance of v2 in the grouped distributions. 

In a similar way one can obtain such formulae for seminvariants as may be 
required. Through the sixth, the formulae for the Sheppard’s corrections for 
the seminvariants calculated from a grouped distribution of a discrete variate are: 


© 1\ 
b+ in — (1 -a)5 


— ls + 3a + os 


. . . 1\ k 

= l, + 3Vo:75 + 431 + 6ro2 + ro + ¢ ore 4) a 
= i+ 10r 11:79 93 + 5a + 10As2 + 10X23 + dos 
le + LV 11:95,» + 10v2:,, a? 3073:1, sie 90r2:+, 20 


ke 
m® } 252° 


= = - ~ = 1 
+ 6As1 + 15a + 20A33 + 15% + ros — ¢ ae ) 


In these formulae, v;;:,,,,, is the 7jth central product moment of v, and », in the 
grouped distributions. 

To illustrate these formulae numerically and to facilitate comparison with 
Pierce’s results, we will use the example he chose. His ungrouped distribution 
was: 
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With origin at vy = 4, we have the following table of moment characteristics 
of these four distributions: 


Distribution| vy; 








9819 | 238849317 | 50388966 | 
60? | 604 


247004154 54 | 


60? 603 608 


10179 567162 | 557840277 | 


$820 1317600 | 528282000 
602 603 | 604 


~ 294904800 





602 603 


| bi 
9606 622440 | 441657198 | 1 ati 
— 
N 
po 


B2 = eo Ms = Az! Ms 


Original | 460 642400 305034000 | 138079200 | 


Distribution’ ~~ 60 603 604 | 604 








From the table, 


— 9606 

602 

x — 622440 
” 608 


- ex? 441657198 
= do + 3r20 = — 60 ss 


We further compute: 


(obvi) _ 6780 _ 
3 603 ” 
_ 8705412 _ 
604 
_ 96774 _ _ _ 2360946 
G04 | 604 
—72978 
604 


Yoo — V2 = 


— 96774 


as on 
vos 3702 = Ere 
) 
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Zva _ 52, _ 330048 
604 


= — — pb = 
3 


_ 163839996 


ls Aso — 3V2:1. = 49 — 39% — Sve. = 60: 
( 1) Ka? 
m/12 3 


1 3\ k' 
(1 =)(7 - 3) sto =? 


With these values one may check the formulae (12) and (13) as far as weight 
four. For example: 


— 9006 , 24 _ 2 10 
ee 6 " 6 3 60? 


M ant (163839996 + 991494 — 34821648 — 437868 — 96774 + 8640000) 


_ 138079200 
60* 


It may appear at first glance that since 
Dre = E[v,(5r;)"] 


and could be expressed by means of the notation, vu:,,,,,‘, the notation in (12) 
and (13) could be made more uniform. It could be but at the expense of greater 
complexity in these two sets of results. Moreover, it is convenient that rs 
is expressible in terms of ».’s in precisely the same way that product semin- 
variants are ordinarily expressible in terms of product moments. 

Pierce’s results differ from the above not only in their mode of derivation but 
also in the fact that they express 7,9’s and l,’s in terms of the characteristics of 
the ungrouped distribution and moments and seminvariants of moments in the 
grouped distributions. Thus as they stand they are not formulae for Sheppard’s 
corrections. 

Finally it must be remarked that in comparison with the usual formulae for 
Sheppard’s corrections, the formulae (10) and (13) introduce quantities the 
magnitudes of which are not known in general except that ordinarily they are 
quite small. It is hoped that results on this point will be forthcoming soon. 
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ON THE ANALYSIS OF VARIANCE IN CASE OF MULTIPLE 
CLASSIFICATIONS WITH UNEQUAL CLASS FREQUENCIES 


By ABRAHAM WALD’ 


Columbia University 


In a previous paper’ the author considered the case of a single criterion of 
classification with unequal class frequencies and derived confidence limits for 
o”/o° where o” denotes the variance associated with the classification, and o° 
denotes the residual variance. The scope of the present paper is to extend 
those results to the case of multiple classifications with unequal class frequencies. 

For the sake of simplicity of notations we will derive the required confidence 
limits in the case of a two-way classification, the extension to multiple classifica- 
tions being obvious. 

Consider a two-way classification with p rows and g columns. Let y be the 
observed variable, and let n;; be the number of observations in the 7th row and 
jth column. Denote by y$} the kth observation on y in the ith row and jth 
column (k = 1, ---,%:;). Let the total number of observations be N. We 
order the N observations and let y. be the ath observation on y in that order. 
Consider the variables: 


t, ti, +++ bp, M1, +++ Uy, 





and denote by ¢. the ath observation on ¢, by tia the ath observation on ¢; and 
by v;2 the ath observation on v;. The values of t., tia and v;. are defined as 
follows: 


tz =1(a=1,---,N), 

tia = 1 if yq lies in the 7th row, 

tin = 0 if y. does not lie in the 7th row, 
Via = 1 if yz. lies in the jth column, 


Vja = 0 if y. does not lie in the jth column. 


We make the assumptions 
(k) (x) 
Wig =U +raty, 
. k ” 2a i 

where the variates a$;, e«,2a,@=1,---,pjg=i,---,qk = 1,--- , m3) 

° ° . : kb) « 2 ° 
are independently and normally distributed, the variance of 2} is o”, the vari- 
ance of ¢ iso”, the variance of 7; is o’”, and the mean values of ¢ and 7; are 
zero. 


1 Research under a grant-in-aid from the Carnegie Corporation of New York. 
2‘‘A note on the analysis of variance with unequal class frequencies,’’ Annals of Math. 
Stat., Vol. 11 (1940). 
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Let the sample regression of y on ¢, t), +++ , tpt, U1, °** 5 Vgea be 


Y — at + bit; + en + bp-ttpr + div. + pliacis + dy 1 . 


We want to derive confidence limits for 


o”/a' =n. 
Let us introduce the notations: 

Du tatia = oi (¢@=1,---,p—1), 
¥ bitin = Gest Gj=1,+--,q¢—D), 
> tiatia = Mi; (i,j = 1, p- t), 
2, itia = Cip143 G@=1,---,p—1;j =1,---,¢—-1)), 
z. Via ja = Ap-14i p—1+i (i,j = 1, »gq— 1), 

II ees |] = | aes [7 (i,j = 0,1,---,p+q-— 2). 


Let the regression of ai on t, th, ++ ,tp-1, U1, °°* , Vg-1 De 

X = att + bith +--+ + dpatpa + din + +++ + dgiave. 
The regression of ¢; + 7; on the same independent variables is evidently equal to 
@li +--+ + eptp + mri + e+ + nada 
= (ng + ep)t + (a — eli t+» + (Qu — €p)tpa 
+ (m — 1g) + e+ + (mga — 
-—tpiandv, =t—v —--+ — v1. Hence 


(1) b; = b; + (& — &), ((=1,---,p— 1), 









Na)Ve-1 5 








sncetp = t—t—-- 





and therefore 





Gd; = Tdzd; + OT (e;—ep) (€j—ep) = Cio" + Fexe; - Fenen 
[cs + (i + 5:;)d' Jo’, (i,j =1,---,p- 1), 
where 6;; is the Kronecker delta, i.e. 5;; = 0 for? ¥ j and 6; = 1. Denote 


G; + (1 + 4;;)\” by c;;. Since the expected value of b; is equal to zero, on 
account of (1) also the expected value of b; is equal to zero. Let 


(2) 


| ges || _ Lees II, (i,j = 1, +++, p— 1). 
Then 
1 p—l p—1 
(3) = p zs gi; bib; 
O° j=l i=l 
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has the x’-distribution with p — 1 degrees of freedom. The expression 
i< 
(4) o dX (Ya me TJ 


has the x’-distribution N — p — q + 1 degrees of freedom. The expressions 
(3) and (4) are independently distributed. Hence 
(5) N-p-qt+l LZgiz bib; 


yo 1 Z(Ya ae Y.)*’ 
has the F-distribution (analysis of variance distribution). We will now show 


that (5) is a monotonic function of \*. It is known that =2g;;bib; is invariant 
under linear transformations, i.e. 








EDgisbib; = VDgs jib; , 


where b; is an arbitrary linear function, say wiabi + +++ + bip-idp-s Of by, «+> 
boi (@ = 1,---,p — 1) and 


’ 


I} gia |] = || ooso5 |. 
We can choose the matrix || y;; || such that 
€; = paler — €p) + +++ + wip-rlepi — €), (Q=1,---,p- 1), 


are independently distributed and o; = o”. The coefficients u;; of course do 
not depend on o’. We have 


ovin; = art's} + diyo", (6;; = Kronecker delta). 


Now let 


bs = vabi + coe + ripedy-1, (vy = 1,---,p— 1), 


° ° ° * *i 
where || »;; || is an orthogonal matrix and is chosen such that b;”, --+ ,bp1 


are independently distributed. On account of the orthogonality of || :; || we 
obviously have 


2 2 2 ° . 
Od; = oD” 4- o’ : O4'0" = 0 for 2 gt je 
Hence 


P p—l1 — 
(6) > LD ggbid; = Do 


in “3 23° 
i=l opr” +rX\o 


The right hand side of (6) is evidently a monotonic function of \” which proves 
our statement. The endpoints of the confidence interval for \” are the roots 
in \” of the equations 





(7) N—p—gt1 Szgibib} _ pp. Na=pr-~aqt1 22gibidi _ p 
p—1 Sye-Yae |’ p—1 Z(ye — Ya)? 


where F, denotes the upper, and F, the lower critical value of F’. 
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The derivation of the required confidence limits in case of classifications in 
more than two ways can be carried out in the same way and I shall merely state 
here the results. 

Consider r criterions of classifications and denote by p, the number of classes 
in the uth classification (u = 1,---,7r). Denote by ni,...;, the number of 
observations which belong to the 7,th class of the first classification, 72th class 
of the second classification, --- , and to the 7,th class of the rth classification. 
Let yf... 4, be the kth observation on y in the set of observations belonging to 
the classes mentioned above (k = 1, --- , ni,-..:,). We make the assumption 

Yor one = Tyrone, H by tooo + a? 
where the variates 


(1) ef” aa ° a . —_ 
tig» €ip 9 9 °° » Gi, (tu = 1 +++, Du; U o 1, ove 9%, k= 1, Niysrig)y 


are independently and normally distributed, the variance of af" ..i, is o°, the 
variance of «“) is o and the mean value of ¢§“’ is zero (i, = 1,+-+, pu; 
= 1, - . 9 ¥). 

Let N the total number of observations. We order the observations in a 
certain order and denote by y, the ath observation in that order (a = 1, --- , N). 
Consider the variables: 


ec”, (u=1,+++ 73 ty = 1, ++, Dus 


and denote by ¢, the ath observation on ¢ and by ¢{“’a the ath observation on 
t{*’. The values of ¢. and ¢{“2 are given as follows: 


= 1 (a = 1, ++ , NV), 


= 1 if y, lies in the 7,th class of the uth classification, 


ti“ = 0 if yq does not lie in the i,th class of the uth classification. 
Let the sample regression of y on t, t{“’ be given by 


=at+ > z bf’ °° 


u=xl iy=l 


tulu 


that 0: = og = --- = 6, =0. The matrix || C$“) || (Qu, j. = 1,-++, pu — 1) 
can be calculated by known methods of the theory of least squares. Let 


Let the covariance of b{“’ and b§“ be given by C{“) o° under the a 


Hossa] = |] CG. (LA Biyiu) Xe [7 (iuy ju = 1, +++, Pu — 2), 
where 6;,;, is the Kronecker delta and \7, = o1/o. Then the lower and upper 
confidence limits for \% are given by the roots in \%, of the equations 


Pu—l pu-l 
N- = Putr—1 dD D giro ose 
(8) se anenrenania a = F, (i = 1, 2), 
” > (Ya — Ya)? 
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where fF is the upper and F;, the lower critical value of the analysis of variance 
distribution with p, — 1 and N — >> p, + r — 1 degrees of freedom. In 
u=l 


case of a single criterion of classification the confidence limits (8) are identical 
with those given in my previous paper. 


THE FREQUENCY DISTRIBUTION OF A GENERAL MATCHING 
PROBLEM 


By T. N. E. GREvILLE 


Bureau of the Census 


1. Introduction. This paper considers the matching of two decks of cards of 
arbitrary composition, and the complete frequency distribution of correct 
matchings is obtained, thus solving a prokhlem proposed by Stevens.’ It is also 
shown that the results can be interpreted in terms of a contingency table. 

Generalizing a problem considered by Greenwood,’ let us consider the matching 
of two decks of cards consisting of ¢ distinct kinds, all the cards of each kind being 
identical. The first or “call” deck will be composed of 2; cards of the first kind, 
72 of the second, etc., such that 


itt tisteee $= 7; 


and the second or “target”? deck will contain 7; cards of the first kind, jo of the 
second, etc., such that 





















AtRptess the Hn 
Any of the 2’s or j’s may be zero. It is desired to calculate, for a given arrange- 
ment of the ‘call’ deck, the number of possible arrangements of the “target” 
deck which will produce exactly r matchings between them (r = 0, 1, 2, --- , n). 
It is clear that these frequencies are independent of the arrangement of the call 
deck. For convenience the call deck may be thought of as arranged so that all 


the cards of the first kind come first, followed by all those of the second kind, 
and so on. 


2. Formulae for the frequencies. Let us consider the number of arrange- 
ments of the target deck which will match the cards in the kith, keth, --- , ksth 
positions in the call deck, regardless of whether or not matchings occur elsewhere. 
Let the cards in these s positions in the call deck consist of c, of the first kind, 
C2 of the second, ete. Then: 


Cte tees $e = 8. 


1 







The number of such arrangements of the target deck is 
(n — s)! 
(1) : peace 
I] (jn — cn)! 


1W.L. Stevens, Annals of Eugenics, Vol. 8 (1937), pp. 238-244. 
2 J. A. GREENWOOD, Annals of Math. Stat., Vol. 9 (1938), pp. 56-59. 
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lor fixed values of the c’s, the s specified positions may be selected in 


t 2 ' 
(2) lieaaitesaiiiie 
h=l1 cn! ‘1. = Cr)! 
ways. 
Consider now the expression 


t 
(n — s)! [J i;! 
(3) V, ssgutninsneemneisenie eines 
TT] cx!(in — cn)! — on)! 
h=1 
obtained by summing the product of (1) and (2) over all sets of values of the 
numbers ¢;, C2 +--+ , ¢; satisfying the conditions: 


t 
0OSaSh, Ch S jny and 2d CG. = & 

Let W. denote the number of arrangements of the target deck which result in 
exactly s matchings. Then it is evident that V, exceeds W, , since the former 
includes those arrangements which give more than s matchings, and these, 
moreover, are counted more than once. Consider an arrangement which 
produces wu matchings, where wu > s. Such an arrangement will be counted 
once in V, for every set of s matchings which can be selected from the total of 
u—that is “C, times. In other words, 


Ve = v, + "HC Waa + ns + ill + "C.Wa. 
It has been shown’ that the solution of these equations is 


(4) W, = V; —_ CT as + as am + (—1)"” "So 


3. Computation of the frequencies. Equations (3) and (4) apparently give 
the solution of the problem, but in practice the labor of carrying out the sum- 
mation indicated in (3) would often be very great. However, (3) may be re- 
written in the form 


(5) =~ 82g. 


t 
Is! 


h=l 


a, iH in! ja! . 


hal Cr! (tn — Cr)! (jn — Cr) 


3H. GEIRINGER, Annals of Math. Stat., Vol. 9 (1938), p. 262. 
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It will be seen that H, is the coefficient of x* in the product 


7 
t th 


in! jn! x" \ 
6) IT ie KV (in — KS Gin — BY? 


t 
where 7, denotes the smaller of 7, and j,. The factor [[ j,! was included in 
h=l 


H, in order to make the coefficients in the polynomials of (6) always integers. 
Equation (4) may now be written in the form 


n —," 
W, = 7 (—1)*” +. o s)! 2. 
vs II is! 


h=el 












(7) wal S(eyTsla-s!, 
‘“ALeE-H =... 
( ) TT in! 


a form which lends itself to actual computation. 


4. Factorial moments. The factorial moments of the frequency distribution 
of the number of matchings are easy to compute. Let m, denote the sth factorial 
moment, so that 


=. r® Ww. 
(8) m, = =. 

>» W,; 
r=0 


Substituting from (4) 


n n n \ 
pi row, _ > i > i (—1)*” "C. V.S 


Tr==8 r=s u=r J 


Reversing the order of summation and simplifying 
a ) 


Lr°w, = {u ¥ 2. t~t" Baal =s!Vi 





n! 


Ve= DW, = - ? 


_ II i! 


h=) 
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and from (5) and (8), 


(10) 


5. Mean and variance. From (6) 
t 
(11) H, = a thyh 
and 


t t 
(12) H, = +2 in(in — 1)jn(jn — 1) + p> integnde « 
hyek 


Hence the mean number of matchings is 


t 
D inj 


n 


(13) 


The variance pe is 


n(n — 1) 


t t 
Mm, + m, — m, = ———— ln dX tn(in — 1)jn(jn — 1) + 2n 2 thir Jnjk 
hek 


t t 2 
+ n(n — 1) dX thjn — (n — (> inn) | ’ 
or 
t t 
(4) m= 1 inin) — tralia +50) + D iis}, 
Gm = D | h=l h=l 
In the special case j: = jo = -+- = Jj: = J, these formulae become 
: t 
. J 2 “2 
M, = = —__ — , 
1 Jy He n(n aa 1) (x dX it) 
These formulae have previously been given by Stevens,‘ and those for the 
special case also by Greenwood. The maximal conditions for the variance, 


given by Greenwood for this particular case, apparently can not be put in a simple 
form for the general case. 


6. Unequal decks. Suppose the call deck contains m cards, m < n, and is to 
be matched with m cards selected from the target deck. It can be assumed 
without loss of generality that the first m cards in any arrangement of the target 
deck are the ones to be used. The formulae of this paper can be applied to this 


‘W.L. Stevens, Annals of Eugenics, loc. cit., Psychol. Review, Vol. 46 (1939), pp. 142-150. 
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more general problem by the expedient of imagining n — m blank cards to be 
added at the end of the call deck and regarding these as an additional kind, 
It is thus apparent that formulae (13) and (14) apply without modification to 
this altered situation. 


7. Application to contingency table. Stevens’ has considered the distribution 
of entries in a contingency table with fixed marginal totals, and has pointed out 
that the problem of matching two decks of cards may be dealt with from that 
standpoint. A contingency table classifies data into n columns and m rows, 
and we may consider the row as indicating the kind of card which occupies 
a given position in the call deck, the columns having the same function with 
respect to the target deck. Stevens defines a quantity c as the sum of entries 
in a prescribed set of cells, subject to the condition that no two cells of the set 
are in the same row or column, and mentions as unsolved the problem of the 
exact sampling distribution of c. 

We now have at our disposal the machinery for solving this problem. Fol- 
lowing Stevens’s notation, let a; , a2, +--+ ,@m denote the fixed row totals and 
b,, be, --+ , b, the fixed column totals, while z,, denotes the frequency of the 


I 
cell in the rth row and the sth column. Then, letc = >> Xr,s, , Where | does 
h=1 


™m n 

not exceed either morn. Imagine two decks of N cards (w =dYia=»> i), 
h=1 h=1 

the first containing a; cards of one kind, a2 of another, etc., and the second 


containing b; cards of one kind, be of another, etc. Moreover, let the r;th kind 
in the first deck and the s,th kind in the second deck be the same kind (h = 
1, 2,--- ,1), the other kinds being all different. Evidently c is the number of 
matchings between the two decks. Hence, the methods of this paper can be 
used to obtain the distribution of c. The formulae we have obtained agree with 
those for the expected value and variance of c given by Stevens. 


ON METHODS OF SOLVING NORMAL EQUATIONS 
By Paut G. Hoe. 


University of California, Los Angeles 


There seems to be considerable disagreement concerning what is the most 
satisfactory method of solving a set of normal equations. Since such informa- 
tion as errors of estimate and significance of results is usually desired in addition 
to the solution, in its broader aspects the problem is one of deciding what is the 
most satisfactory method of calculating the inverse of a symmetric matrix. 

For equations with several unknowns some compact systematic method of 


’W. L. Stevens, Annals of Eugenics, loc. cit. 
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calculation is necessary to eliminate much of the labor involved in the ordinary 
method of calculating the inverse from its definition. Among the more common 
of such systematic methods are those associated with the names of Chio,’ Gauss,’ 
Doolittle,’ and Aitken.’ In addition, A. A. Albert* recently called attention 
to a method implicit in elementary matrix theory. There are also various 
iterative schemes, and schemes which are but slight variations of the above 
methods. In this note only the methods associated with the above names will 
be considered, and for convenience they will be labeled with those names, regard- 
less of who should be given credit for them. 

The purpose of this note is to show that when the calculation of the inverse is 
systematized, all of the above methods are fundamentally equivalent and merely 
involve a different arrangement of work. Consequently, any advantage in calcu- 
lating time for any particular method will arise through such features as a 
simpler technique or less copying, rather than through fewer multiplications and 
divisions. 

By the method of Chio is meant the evaluation of determinants by the pivotal 
method of reduction. Since all of the methods mentioned above use pivotal 
reduction, the method of Chio will not be treated as a distinct method. Fur- 
thermore, since Gauss’ method is incorporated in that of Aitken, it will be neces- 
sarv to consider only the methods of Aitken, Doolittle, and Albert as distinct. 

First consider the method of Albert, which is based on the following matrix 
properties. Let the matrix A be subjected to a sequence of row transformations 
leading to the matrix A’. Then, writing A = IA, it follows from a theorem in 
matrix theory that A’ = I’A, and consequently that A‘A' = I’. If row trans- 
formations are chosen which make A’ = I, then A’ = I’. This states 
that if the same row transformations are applied to the identity matrix as were 
used to reduce A to the identity matrix, then the resulting matrix will be the 
desired inverse. The customary manner of reducing A to I is to work for zeros 
in columns as follows: 




















an 1 a2 Qin | 
1) an ay 
— Qi2 *** Gin | 
| a21 a21 | 
|}@21 ez +++ Gen 0 ae — — “77 1” ade 
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1 | 


1See, for example, Whittaker and Robinson, The Calculus of Observations, p. 71 and p. 
234. 

*See, for example, Croxton and Cowden, Applied General Statistics, 1939, p. 716. 

3 Roy. Soc. Edin. Proc., Vol. 57 (1936-37), p. 172. 

4Am. Math. Monthly, Vol. 48, No. 3 (1941), p. 198. 
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where new letters are introduced for new elements after each reduction. After 
zeros are obtained below the main diagonal, zeros are obtained above the 
diagonal by starting with the last column. If now these operations are per- 
formed in the same order on I, the result will be A. 

Next consider the method of Aitken, which is based on the evaluation of a 
bordered determinant, namely, 
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To obtain A it is merely necessary to evaluate determinants of this type and 
divide them by | A |. Aitken’s method evaluates all such determinants simulta- 
neously, using Chio’s reduction technique in much the same manner as illustrated 
above with Albert’s method. Thus, 
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When zeros are obtained below the main diagonal to the left of the vertical 
dividing line, the matrix in the lower right section will be A’. This follows from 
the fact that the elements of this matrix will be the evaluations of bordered 
determinants, like those of the previous paragraph, divided by ane --- = | A |. 

It will be observed that the operations on A in Albert’s method which produce 
zeros below the main diagonal are the same as those which occur above the hori- 
zontal dividing line in Aitken’s method. This set of operations is performed 
simultaneously on I, since the upper right section of Aitken’s scheme is I. Fur- 
thermore, obtaining a zero for an element below the horizontal line and to the 
left of the vertical line, is equivalent to obtaining a zero for the clement corre- 


bo 
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sponding to the same row and column in the section above the horizontal, pro- 
vided the preeeding columns contain zeros above the diagonal. But obtaining 
zercs above the main diagonal of A constitutes the second set of operations in 
Albert’s method to obtain A’ = I. Thus, the operations in Aitken’s method 
which produce zeros in a given column for elements above the horizontal line 
are merely the first set of operations in Albert’s method, while those which 
produce zeros below the horizontal line are the second set of operaticns in reverse 
order. Sinee, in Aitken’s scheme, the first set of operations is performed on I 
in the upper right section and the results are transferred a row at a time to the 
lower right section, where they are in turn operated upon by the second set of 
operations, this lower right section is merely I operated upon by the entire set 
of operations of Albert’s method. Consequently, Aitken’s and Albert’s methods 
are the same except for the order in which operations are performed and differ- 
ences arising therefrom. Since Aitken’s method performs these operations more 
compactly, it is to be preferred to that of Albert. 

Next consider the method of Doolittle, which is described by following 
the instructions given in the first column in the table shown on page 348. 
The forward solution is completed after such sectional operations. For a 
given k column, the backward solution is obtained as usual by substitution in 
the last row of each section taken in reverse order. 

If all summations in each section are performed in pairs and the sums recorded 
each time, rather than being performed in one operation, the forward solution 
of the Doolittle method will be found to be a rearrangement of the work occurring 
above the horizontal line in Aitken’s method. Thus the first lines of each 
section give the matrix above the horizontal line in Aitken’s scheme. Then, 
except for signs, J’ and the sums of the first two lines of the remaining sections 
give the result of Aitken’s first sequence of operations above the horizontal. 
Then, except for signs, JJ’ and the sums of the first three lines of the remaining 
sections give the result of Aitken’s second sequence of operations above the 
horizontal, ete. 

The back solution involves precisely the same operations as those making up 
the second set of Albert’s sequence of operations to obtain zeros above the main 
diagonal. Since these were shown to be a rearrangement of operations in 
Aitken’s method, it follows that the methods of Aitken and Doolittle are the 
same except for the order of operations and differences arising therefrom. Hence 
all three methods are basically the same when systematized for a calculating 
machine. 

Because of this equivalence, the number of necessary multiplications and 
divisions will be the same for all three methods, and will be found to be 
1n?(n + 1). Since Aitken’s method is to be preferred to that of Albert, it will 
suffice to compare the methods of Aitken and Doolittle for calculating con- 
venience. 

The Doolittle method possesses several distinct advantages. First, its multi- 
plications occur a row at a time with one of the factors constant for that row; 
consequently the keyboard remains unchanged for a given row of operations. 
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Aitken’s method, however, consists of calculating successive cross products’ 
which requires clearing of the keyboard after each such operation. Secondly, 
there are fewer additions in the Doolittle method. It sums 7 quantities at a 
time in section 2, while Aitken’s cross products always involve the sum of two 
quantities. Because of the necessity of calculating the complements of negative 
sums, this difference becomes important when the number of variables is large. 
A third feature in favor of the Doolittle method is the ease of performing the 
calculations without previous experience. It may be easier to understand how 
to calculate cross products, but actually the calculations of the Doolittle method 
are easier to perform. Aitken’s method requires some experience with it, if one 
is to avoid repeating certain calculations which would result from calculating all 
cross products mechanically. The comparative amount of copying in the two 
methods depends upon the number of variables involved. 

From the above considerations, it may be concluded that the Doolittle method 
is to be preferred among those considered in this paper for solving a set of normal 
equations or calculating the inverse of a symmetric matrix. However, if a 
single calculating technique is desired which can be used for nonsymmetrical 
equations as well, then the method of Aitken is to be preferred. 
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CONDITIONS THAT THE ROOTS OF A POLYNOMIAL BE LESS THAN 
UNITY IN ABSOLUTE VALUE 


By Pau A. SAMUELSON 
Massachusetis Institute of Technology 


1. Introduction. In econometric business cycle analysis, probability the- 
ory, and numerical mathematical computation the problem of convergence of 
repeated iterations arises. The solution of the difference equations defining 
such a process can in a wide variety of cases be shown to be stable in the sense of 
converging to a limit if a certain associated polynomial 


(1) f(x). = pox” + pit” + --» + pn = 0, 


has roots whose moduli are all less than unity. 


Thus, for ‘‘timeless” linear difference equation systems of the most general 
type, convertible into normal form, 


(2) Q(t+1) = . ai;Q;(t), 


the polynomial is the characteristic or determinantal equation, 
(3) f(z) = | ai; — 28s; | = 0, 


which when expanded out is of the form (1). The roots of this equation, when 


multiplied by suitable polynomials in ¢, give the exact solution of the problem 
in the form 


(4) Q(t) = 3 gi(t)axi , 


where m is the number of distinct roots, and the g’s are polynomials of degree 
one less than the multiplicity of the respective root. If complex roots occur, 
they do so in conjugate pairs and can be combined to form damped, undamped, 
or anti-damped harmonic terms. All terms go to zero as t approaches infinity 
if, and only if, the absolute value of each z is less than unity. 

For non-linear systems the exact solution does not take this form, but in the 
neighborhood of an equilibrium point the roots of an associated polynomial, 
except in singular cases, do determine the stability of the system. 

As far as the writer is aware, there does not appear in the literature an account 
of necessary and sufficient conditions for the roots of a polynomial to be less than 
unity in absolute value. This is in contrast to a related problem which arises 
in connection with the investigation of stability of dynamical systems defined by 
differential equations. These have associated with them a polynomial whose 
roots provide solutions in the form 


(5) gi(t)e** ; 
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or for non-linear systems infinite power series in such terms. It is required, 
therefore, to determine complete conditions under which the real parts of all 
roots must be negative. ° 

This problem has been solved by Routh’ in a manner which leaves little to be 
desired. Determinantal expression of his conditions in a slightly modified form 
was made by Hurwitz’ who apparently was unaware of Routh’s work, and by 
Frazer and Duncan* who were unaware of the Hurwitz results. A brief outline 
of Routh’s mode of attack will prove instructive in dealing with the problem 
at hand. 





2. Routhian analysis of sign of real parts of roots. Routh realized 
that the condition that all coefficients be positive—the leading coefficient having 
been made so—was necessary, but not sufficient unless all the roots were real. 
But a “derived” equation of degree n(n — 1)/2 whose roots equal the sums of the 
roots of the origina] equation taken two at a time has real roots which are simple 
sums of the real parts of those of the original equation. In consequence, it is 
necessary and sufficient that the coefficients of the original and the ‘‘derived”’ 
equation all be positive. 

Thus, valid necessary and sufficient conditions are presented. However, they 
are disadvantageous from two points of view. First, they are not all independ- 
ent, being n(n + 1)/2 conditions in number, whereas only n are necessary. Sec- 
ondly, despite several ingenious methods devised by Routh, it is not easy to 
compute them in the general case. 

Recognizing these difficulties, he therefore began anew from an entirely 
different angle. Utilizing a theorem of Cauchy concerning the relationship 
between the behavior of a polynomial on a closed contour in the complex domain 
and the number of roots within that closed curve, he derived necessary and suffi- 
cient conditions, which may be written in the slightly more convenient deter- 
minantal form of Hurwitz and Frazer and Duncan as follows: 


Pi Ds! 


Po Pe 
P3-** Pos-1 


To = po > 0; T, = pi > 0, T, = > 0, 




















+ P2s—2 


1}. J. Routh, A Treatise on the Stability of a Given State of Motion, (London, 1877), 
Chaps. 2 and 3; Advanced Rigid Dynamics, 6th ed., London, 1905, Chap. 6. 

? Hurwitz, Math. Ann., Vol. 46 (1895), p. 521. 

*R. A. Frazer and W. J. Duncan, Royal Soc. Proc., Series A, Vol. 124 (1929), p. 642. 
Also R. A. Frazer, W. J. Dunean. and A. R. Collar, Elementary Matrices, Cambridge Uni- 
versity Press, 1938, pp. 151-155. 
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The law of formation of these determinants is obvious. In the first row the 
odd p’s starting with the first are listed. Within each column the p’s diminish 
one unit at a time. Any p with negative subscript derived by this formula is 
treated as zero, and all p’s of subscript higher than the degree of the equation 
are set equal to zero. With this convention, for p) made positive, complete and 
independent necessary conditions are that all principal minors of 7’, formed by 
deleting successively the last row and column must be positive. These condi- 
tions are n in number and are independent. 


3. Complete, independent, necessary and sufficient conditions. Corre- 
sponding to Routh’s first attack on the problem, we might consider an equation 
of degree n(n — 1)/2 whose roots equal the products two at a time of the original 
equation’s. If this equation and the original equation have real roots less than 
unity in absolute value, our problem is solved. This is guaranteed if, and only 
if, two further transformed equations with roots equal to the squares minus unity 
of the roots of the original and derived equations respectively all have positive 
coefficients. These conditions are necessary and sufficient, but not independent, 
and cannot be easily computed in the general case. Therefore, I follow Routh’s 
example and approach the problem from a different point of view. 

When the roots of f(x) = 0 are plotted in the complex plane, they must all lie 
within the unit circle if their absolute values are to be less than unity, and con- 
versely. We might therefore attempt to apply Cauchy’s theorem. However, 
it is not necessary todo so. Routh has shown what the conditions are that there 
be no roots in the right-hand half-plane. Can we find a complex transformation 
of variables which carries the unit circle into the left-hand half-plane? 

The answer is in the affirmative. The linear complex transformation 


(7) ne? oe ees 


z— 1’ z—1 


will accomplish this. But after substituting for x its value in terms of z, we 
cease to have a polynomial but rather a rational function of z as follows: 


Y pile +1)" — 1) 
(8) fla) = 5(7 : D Al cccmaecisiesesiicnsns stl 


z—1 (2 — 1)" 


We need only consider the polynomial in the numerator, i.e., 
(9) g(z) = X m2” =0. 


In order that the roots of the original equation be less than unity, in absolute value, 
it is necessary and sufficient that the real parts of the roots of equation (9) be negative. 
Once we determine the coefficients (7;) in terms of the original p’s, we can easily 
apply Routh’s theorems. This yields n + 1 necessary and sufficient conditions, 
all of which are independent. 
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Expanding the numerator of the right-hand side of (8) and collecting terms, 
the following explicit formulas for the ’s are directly obtained: 


(10) = a. Pi y n—jCi-~(—1)* iCi, 


7 =0 k= 


where 


v! 


ve 


io — wie’ 
and 
m(i, 7) = the smaller of 7 and j. 


For fourth and higher degree equations literal substitution, while always 
possible, results in complicated expressions. It is preferable, therefore, to com- 
pute the w’s numerically and then apply the conditions of (6) directly. 

Other necessary conditions can be easily derived, but they will be dependent 
upon these. Thus, each 7 must be positive; but this is not, by itself, sufficient. 
Or, adding 7» and z, we find 


(11) T + tT, = Po t+ Po +mwat--- > O, 


i.e., the sum of the even p’s must be positive. Similarly, still other linear sums 
of other 7’s will result in cancellation of certain of the p’s. Except on special 
occasions there is probably no labor saved by utilizing conditions derived in 
this way. 

One obvious but useful necessary condition will be stated without proof. 
If one forms polynomials from subsets of the coefficients of a given “stable’’ 
polynomial formed by arbitrary “cuts” which leave adjacent coefficients in 
unchanged order and introduce no gaps within each set, then the resulting poly- 
nomials will all be stable. 

Special sufficiency conditions also can be developed. Carmichael‘ presents 
certain inequalities between the absolute values of the largest root and the coeffi- 
cients of the original equation. For special problems these may be fruitfully 
applied. 


4. Example. In conclusion I apply the conditions derived here to a well- 
known numerical equation determined statistically by Tinbergen’ in the analysis 
of economic fluctuations. It is a fourth order difference equation with constant 
coefficients, 


(12) Zi — .898Z141 + .220Z:-2 — .013Z:-3 — .027Z.4 = 0 


*R. D. Carmichael, Amer. Math. Soc. Bull., Vol. 24 (1918), pp. 286-296. 


5 J. Tinbergen, Business Cycles in the United States, 1919-1982, League of Nations, 1939, 
p. 140. 
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with the associated indicial equation 
(13) f(z) = x* — .3982° + .22027 — .013x — .027 = 0. 


Its roots have been computed and are known to be less than unity in absolute 
value. This may be verified by computing 


0:782 > 0 





™ = 3.338 > 0 
™ = 5.398 > 0 
(14) ™ = 4.878 > 0 
4 = 1.604 > 0 


T2 = 14.204 > 0 
T; = 43.177 > 0 


To compute the same results by cross-multiplication the work is arranged as 
follows: 






Wo 
.782 
m1 
(15) 3.338 4.878 
W182 — Worms 1374 — 0 
14.204 7.824 
3(%1%2 — Wows) — Wiw3%4 


43.177 

























It may be remarked that the presence of a negative coefficient anywhere in 
the table is an immediate indication of instability, and that there is no necessity 
to continue the computation until a negative sign appears in a leading coefficient. 
This fact often saves much labor. 


VALUES OF MILLS’ RATIO OF AREA TO BOUNDING ORDINATE AND OF 
THE NORMAL PROBABILITY INTEGRAL FOR LARGE VALUES 
OF THE ARGUMENT 


By Rosert D. Gordon 






Scripps Institution of Oceanography 


A pair of simple inequalities is proved which constitute upper and lower 
bounds for the ratio R,’, valid for zx > 0. The writer has failed to encounter 
these inequalities in the literature, hence it seems worthwhile to present them 
for whatever value they may have. 





1J. P. Mills, ‘Table of ratio: area to bounding ordinate, for any portion of the normal 
curve.” Biometrika Vol. 18 (1926) pp. 395-400. Also Pearson’s tables, Part II, Table IIT, 



















MILLS’ RATIO 


The funetion R, is defined by 


(1) , a | co" é. 


The following relations between R = R, and its derivatives are easily established 
by direct differentiations and substitutions: 


sare vent Ste, ?, 


dx? z dx 


aR d a. 
dx =(1 tay i)* dx? 


Also by ordinary rules 
R, > 0, 


lim «Rk, = 
1°. Suppose that at any point z,; > 0, 1.R > 1. Then by (2) dR/dzx > 0, 
and R, would continue to increase with increasing x: still more, xR, would con- 
tinue to increase, hence we should have zR, > 1 for x 2 2, which contradicts 
(6). Therefore we find eR, S 1 for x > 0, and 


which establishes the required upper inequality. 

2°. Suppose that at any point x. > 0, d’R/dx? <0. Then by (4) d@R/dz* = 
(d/dx)(d°R/dr) < 0 at this point. Since these derivatives are continuous this 
implies that for all x > a2, dR/dx” < [d°R/dz’|,-2, < 0. Then we get the 
inequalities, for 2 > ze 


UR dR [aR dR 

a” E | ee 2) | oe ” [I 

R < R,, + (x — 2) | + 3(x — a)” Fa 
lx |e dx? |e 


where [ ]> indicates evaluation at x = x2. Since [d’R/dz’}, < 0, this implies 
that for sufficiently large x, R, < 0, which contradicts (5). It follows thep 
that (3) is positive, and substitution of (2) gives 


(8) iat 


~ gtd’ 
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We combine (7) and (8) in the double inequality: 


. gp¢* 


9 es ’ 
(9) z+ 1 ~~ B 


This gives for the probability integral the corresponding inequality 


XL 1 —.2/9 1 —t2/2 1 1 —z2/2 
(10) a asf CO eg ts nc. 
r+ i V 27 V/2r x x V/ 20 


It can easily be shown (for 2 > 0) that equalities in (9) and (10) are impossible. 








